(git:ccc2433)
xc_xpbe_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 Calculates the exchange energy for the pbe hole model in a truncated
10 !> coulomb potential, considering the long range part only. Can be used
11 !> as longrange correction to a truncated Hartree Fock calculation
12 !> \par History
13 !> Manuel Guidon (01.2009) : initial version
14 !> \author Manuel Guidon (01.2009)
15 ! **************************************************************************************************
16 
17 MODULE xc_xpbe_hole_t_c_lr
18  USE input_section_types, ONLY: section_vals_type,&
19  section_vals_val_get
20  USE kinds, ONLY: dp
21  USE mathconstants, ONLY: euler,&
22  pi,&
23  rootpi
24  USE xc_derivative_desc, ONLY: deriv_norm_drho,&
25  deriv_norm_drhoa,&
26  deriv_norm_drhob,&
27  deriv_rho,&
28  deriv_rhoa,&
29  deriv_rhob
30  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
31  xc_dset_get_derivative
32  USE xc_derivative_types, ONLY: xc_derivative_get,&
33  xc_derivative_type
34  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
35  USE xc_rho_set_types, ONLY: xc_rho_set_get,&
36  xc_rho_set_type
37 #include "../base/base_uses.f90"
38 
39  IMPLICIT NONE
40 
41  PRIVATE
42 
43 ! *** Global parameters ***
44 
45  PUBLIC :: xpbe_hole_t_c_lr_lda_eval, xpbe_hole_t_c_lr_lda_info, &
46  xpbe_hole_t_c_lr_lda_calc_1, xpbe_hole_t_c_lr_lda_calc_2, &
47  xpbe_hole_t_c_lr_lsd_eval, xpbe_hole_t_c_lr_lsd_info
48 
49  REAL(KIND=dp), PARAMETER :: a1 = 0.00979681_dp, &
50  a2 = 0.04108340_dp, &
51  a3 = 0.18744000_dp, &
52  a4 = 0.00120824_dp, &
53  a5 = 0.0347188_dp
54  REAL(KIND=dp), PARAMETER :: a = 1.0161144_dp, &
55  b = -0.37170836_dp, &
56  c = -0.077215461_dp, &
57  d = 0.57786348_dp, &
58  e = -0.051955731_dp, &
59  f2 = 0.47965830_dp, &
60  f1 = 6.4753871_dp, &
61  clda = -0.73855876638202240588423_dp
62  REAL(KIND=dp), PARAMETER :: expcutoff = 700.0_dp
63  REAL(KIND=dp), PARAMETER :: smax = 8.572844_dp, &
64  sconst = 18.79622316_dp, &
65  scutoff = 8.3_dp
66  REAL(KIND=dp), PARAMETER :: gcutoff = 0.08_dp, &
67  g1 = -0.02628417880_dp/e, &
68  g2 = -0.07117647788_dp/e, &
69  g3 = 0.08534541323_dp/e, &
70  g4 = 0.0_dp
71  REAL(KIND=dp), PARAMETER :: wcutoff = 14.0_dp
72 
73  REAL(KIND=dp), PARAMETER :: eps_rho = 0.00000001_dp
74 
75  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xpbe_hole_t_c_lr'
76 
77 CONTAINS
78 
79 ! **************************************************************************************************
80 !> \brief returns various information on the functional
81 !> \param reference string with the reference of the actual functional
82 !> \param shortform string with the shortform of the functional name
83 !> \param needs the components needed by this functional are set to
84 !> true (does not set the unneeded components to false)
85 !> \param max_deriv controls the number of derivatives
86 !> \par History
87 !> 01.2009 created [mguidon]
88 !> \author mguidon
89 ! **************************************************************************************************
90  SUBROUTINE xpbe_hole_t_c_lr_lda_info(reference, shortform, needs, max_deriv)
91  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
92  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
93  INTEGER, INTENT(out), OPTIONAL :: max_deriv
94 
95  IF (PRESENT(reference)) THEN
96  reference = "{LDA version}"
97  END IF
98  IF (PRESENT(shortform)) THEN
99  shortform = "{LDA}"
100  END IF
101  IF (PRESENT(needs)) THEN
102  needs%rho = .true.
103  needs%norm_drho = .true.
104  END IF
105  IF (PRESENT(max_deriv)) max_deriv = 1
106 
107  END SUBROUTINE xpbe_hole_t_c_lr_lda_info
108 
109 ! **************************************************************************************************
110 !> \brief returns various information on the functional
111 !> \param reference string with the reference of the actual functional
112 !> \param shortform string with the shortform of the functional name
113 !> \param needs the components needed by this functional are set to
114 !> true (does not set the unneeded components to false)
115 !> \param max_deriv controls the number of derivatives
116 !> \par History
117 !> 01.2009 created [mguidon]
118 !> \author mguidon
119 ! **************************************************************************************************
120  SUBROUTINE xpbe_hole_t_c_lr_lsd_info(reference, shortform, needs, max_deriv)
121  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
122  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
123  INTEGER, INTENT(out), OPTIONAL :: max_deriv
124 
125  IF (PRESENT(reference)) THEN
126  reference = "{LSD version}"
127  END IF
128  IF (PRESENT(shortform)) THEN
129  shortform = "{LSD}"
130  END IF
131  IF (PRESENT(needs)) THEN
132  needs%rho_spin = .true.
133  needs%norm_drho_spin = .true.
134  END IF
135  IF (PRESENT(max_deriv)) max_deriv = 1
136 
137  END SUBROUTINE xpbe_hole_t_c_lr_lsd_info
138 
139 ! **************************************************************************************************
140 !> \brief evaluates the pbe-hole exchange in a truncated coulomb potential
141 !> \param rho_set the density where you want to evaluate the functional
142 !> \param deriv_set place where to store the functional derivatives (they are
143 !> added to the derivatives)
144 !> \param order degree of the derivative that should be evaluated,
145 !> if positive all the derivatives up to the given degree are evaluated,
146 !> if negative only the given degree is calculated
147 !> \param params input parameters (scaling,cutoff)
148 !> \par History
149 !> 01.2009 created [Manuel Guidon]
150 !> \author Manuel Guidon
151 !> \note
152 ! **************************************************************************************************
153  SUBROUTINE xpbe_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
154 
155  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
156  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
157  INTEGER, INTENT(IN) :: order
158  TYPE(section_vals_type), POINTER :: params
159 
160  CHARACTER(len=*), PARAMETER :: routinen = 'xpbe_hole_t_c_lr_lda_eval'
161 
162  INTEGER :: handle, npoints
163  INTEGER, DIMENSION(2, 3) :: bo
164  REAL(kind=dp) :: epsilon_rho, r, sx
165  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
166  POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
167  rho
168  TYPE(xc_derivative_type), POINTER :: deriv
169 
170  CALL timeset(routinen, handle)
171 
172  CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
173  CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=r)
174 
175  CALL xc_rho_set_get(rho_set, rho=rho, &
176  norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
177  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
178 
179  dummy => rho
180 
181  e_0 => dummy
182  e_rho => dummy
183  e_ndrho => dummy
184 
185  IF (order >= 0) THEN
186  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
187  allocate_deriv=.true.)
188  CALL xc_derivative_get(deriv, deriv_data=e_0)
189  END IF
190  IF (order >= 1 .OR. order == -1) THEN
191  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
192  allocate_deriv=.true.)
193  CALL xc_derivative_get(deriv, deriv_data=e_rho)
194  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
195  allocate_deriv=.true.)
196  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
197  END IF
198  IF (order > 1 .OR. order < -1) THEN
199  cpabort("derivatives bigger than 1 not implemented")
200  END IF
201 
202  IF (r == 0.0_dp) THEN
203  cpabort("Cutoff_Radius 0.0 not implemented")
204  END IF
205 
206 !$OMP PARALLEL DEFAULT(NONE) &
207 !$OMP SHARED(npoints, order, rho, norm_drho, e_0, e_rho) &
208 !$OMP SHARED(e_ndrho, epsilon_rho) &
209 !$OMP SHARED(sx, r)
210 
211  CALL xpbe_hole_t_c_lr_lda_calc(npoints, order, rho=rho, norm_drho=norm_drho, &
212  e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
213  epsilon_rho=epsilon_rho, &
214  sx=sx, r=r)
215 
216 !$OMP END PARALLEL
217 
218  CALL timestop(handle)
219 
220  END SUBROUTINE xpbe_hole_t_c_lr_lda_eval
221 
222 ! **************************************************************************************************
223 !> \brief intermediate level routine in order to decide which branch has to be
224 !> taken
225 !> \param npoints number of points on the grid
226 !> \param order order of the derivatives
227 !> \param rho value of density on the grid
228 !> \param norm_drho value of gradient on the grid
229 !> \param e_0 derivatives of the energy on the grid
230 !> \param e_rho derivatives of the energy on the grid
231 !> \param e_ndrho derivatives of the energy on the grid
232 !> \param epsilon_rho cutoffs
233 !> \param sx functional parameters
234 !> \param R functional parameters
235 !> \par History
236 !> 01.2009 created [Manuel Guidon]
237 !> \author Manuel Guidon
238 !> \note For numerical reasons there are two different branches
239 !>
240 ! **************************************************************************************************
241  SUBROUTINE xpbe_hole_t_c_lr_lda_calc(npoints, order, rho, norm_drho, e_0, e_rho, e_ndrho, &
242  epsilon_rho, sx, R)
243 
244  INTEGER, INTENT(in) :: npoints, order
245  REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, norm_drho, e_0, e_rho, e_ndrho
246  REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, r
247 
248  INTEGER :: ip
249  REAL(dp) :: my_ndrho, my_rho
250  REAL(kind=dp) :: ss, ss2, sscale, t1, t2, t3, t4, t6, t7, &
251  t8
252 
253 !$OMP DO
254 
255  DO ip = 1, npoints
256  my_rho = max(rho(ip), 0.0_dp)
257  IF (my_rho > epsilon_rho) THEN
258  my_ndrho = max(norm_drho(ip), epsilon(0.0_dp)*1.e4_dp)
259  ! *** Do some precalculation in order to catch the correct branch afterwards
260  sscale = 1.0_dp
261  t1 = pi**2
262  t2 = t1*my_rho
263  t3 = t2**(0.1e1_dp/0.3e1_dp)
264  t4 = 0.1e1_dp/t3
265  t6 = my_ndrho*t4
266  t7 = 0.1e1_dp/my_rho
267  t8 = t7*sscale
268  ss = 0.3466806371753173524216762e0_dp*t6*t8
269  IF (ss > scutoff) THEN
270  ss2 = ss*ss
271  sscale = (smax*ss2 - sconst)/(ss2*ss)
272  END IF
273  IF (ss*sscale > gcutoff) THEN
274  CALL xpbe_hole_t_c_lr_lda_calc_1(e_0(ip), e_rho(ip), e_ndrho(ip), &
275  my_rho, &
276  my_ndrho, sscale, sx, r, order)
277  ELSE
278  CALL xpbe_hole_t_c_lr_lda_calc_2(e_0(ip), e_rho(ip), e_ndrho(ip), &
279  my_rho, &
280  my_ndrho, sscale, sx, r, order)
281  END IF
282  END IF
283  END DO
284 
285 !$OMP END DO
286 
287  END SUBROUTINE xpbe_hole_t_c_lr_lda_calc
288 
289 ! **************************************************************************************************
290 !> \brief evaluates the pbe-hole exchange in a truncated coulomb potential
291 !> \param rho_set the density where you want to evaluate the functional
292 !> \param deriv_set place where to store the functional derivatives (they are
293 !> added to the derivatives)
294 !> \param order degree of the derivative that should be evaluated,
295 !> if positive all the derivatives up to the given degree are evaluated,
296 !> if negative only the given degree is calculated
297 !> \param params input parameters (scaling,cutoff)
298 !> \par History
299 !> 01.2009 created [Manuel Guidon]
300 !> \author Manuel Guidon
301 !> \note
302 ! **************************************************************************************************
303  SUBROUTINE xpbe_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params)
304 
305  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
306  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
307  INTEGER, INTENT(IN) :: order
308  TYPE(section_vals_type), POINTER :: params
309 
310  CHARACTER(len=*), PARAMETER :: routinen = 'xpbe_hole_t_c_lr_lsd_eval'
311 
312  INTEGER :: handle, npoints
313  INTEGER, DIMENSION(2, 3) :: bo
314  REAL(kind=dp) :: epsilon_rho, r, sx
315  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
316  POINTER :: dummy, e_0, e_ndrhoa, e_ndrhob, e_rhoa, &
317  e_rhob, norm_drhoa, norm_drhob, rhoa, &
318  rhob
319  TYPE(xc_derivative_type), POINTER :: deriv
320 
321  CALL timeset(routinen, handle)
322 
323  CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
324  CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=r)
325 
326  CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
327  norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, local_bounds=bo, rho_cutoff=epsilon_rho)
328  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
329 
330  dummy => rhoa
331 
332  e_0 => dummy
333  e_rhoa => dummy
334  e_rhob => dummy
335  e_ndrhoa => dummy
336  e_ndrhob => dummy
337 
338  IF (order >= 0) THEN
339  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
340  allocate_deriv=.true.)
341  CALL xc_derivative_get(deriv, deriv_data=e_0)
342  END IF
343  IF (order >= 1 .OR. order == -1) THEN
344  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
345  allocate_deriv=.true.)
346  CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
347  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
348  allocate_deriv=.true.)
349  CALL xc_derivative_get(deriv, deriv_data=e_rhob)
350  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
351  allocate_deriv=.true.)
352  CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
353  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
354  allocate_deriv=.true.)
355  CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
356  END IF
357  IF (order > 1 .OR. order < -1) THEN
358  cpabort("derivatives bigger than 1 not implemented")
359  END IF
360 
361  IF (r == 0.0_dp) THEN
362  cpabort("Cutoff_Radius 0.0 not implemented")
363  END IF
364 
365 !$OMP PARALLEL DEFAULT(NONE) &
366 !$OMP SHARED(npoints, order, rhoa, norm_drhoa, e_0, e_rhoa) &
367 !$OMP SHARED(epsilon_rho, sx, r) &
368 !$OMP SHARED(rhob, norm_drhob, e_rhob, e_ndrhoa, e_ndrhob)
369 
370  CALL xpbe_hole_t_c_lr_lsd_calc(npoints, order, rho=rhoa, norm_drho=norm_drhoa, &
371  e_0=e_0, e_rho=e_rhoa, e_ndrho=e_ndrhoa, &
372  epsilon_rho=epsilon_rho, sx=sx, r=r)
373  CALL xpbe_hole_t_c_lr_lsd_calc(npoints, order, rho=rhob, norm_drho=norm_drhob, &
374  e_0=e_0, e_rho=e_rhob, e_ndrho=e_ndrhob, &
375  epsilon_rho=epsilon_rho, sx=sx, r=r)
376 
377 !$OMP END PARALLEL
378 
379  CALL timestop(handle)
380 
381  END SUBROUTINE xpbe_hole_t_c_lr_lsd_eval
382 
383 ! **************************************************************************************************
384 !> \brief intermediate level routine in order to decide which branch has to be
385 !> taken
386 !> \param npoints number of points on the grid
387 !> \param order order of the derivatives
388 !> \param rho density on the grid
389 !> \param norm_drho gradient on the grid
390 !> \param e_0 derivatives of the energy on the grid
391 !> \param e_rho derivatives of the energy on the grid
392 !> \param e_ndrho derivatives of the energy on the grid
393 !> \param epsilon_rho cutoffs
394 !> \param sx functional parameters
395 !> \param R functional parameters
396 !> \par History
397 !> 01.2009 created [Manuel Guidon]
398 !> \author Manuel Guidon
399 !> \note For numerical reasons there are two different branches. This code calls
400 !> the lda version applying spin scaling relations
401 ! **************************************************************************************************
402  SUBROUTINE xpbe_hole_t_c_lr_lsd_calc(npoints, order, rho, norm_drho, e_0, e_rho, e_ndrho, &
403  epsilon_rho, sx, R)
404 
405  INTEGER, INTENT(in) :: npoints, order
406  REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, norm_drho, e_0, e_rho, e_ndrho
407  REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, r
408 
409  INTEGER :: ip
410  REAL(dp) :: my_ndrho, my_rho
411  REAL(kind=dp) :: e_tmp, ss, ss2, sscale, t1, t2, t3, t4, &
412  t6, t7, t8
413 
414 !$OMP DO
415 
416  DO ip = 1, npoints
417  my_rho = max(2.0_dp*rho(ip), 0.0_dp)
418  IF (my_rho > epsilon_rho) THEN
419  my_ndrho = max(2.0_dp*norm_drho(ip), 0.0_dp)
420 
421  ! *** Do some precalculation in order to catch the correct branch afterwards
422  sscale = 1.0_dp
423  t1 = pi**2
424  t2 = t1*my_rho
425  t3 = t2**(0.1e1_dp/0.3e1_dp)
426  t4 = 0.1e1_dp/t3
427  t6 = my_ndrho*t4
428  t7 = 0.1e1_dp/my_rho
429  t8 = t7*sscale
430  ss = 0.3466806371753173524216762e0_dp*t6*t8
431  IF (ss > scutoff) THEN
432  ss2 = ss*ss
433  sscale = (smax*ss2 - sconst)/(ss2*ss)
434  END IF
435  e_tmp = 0.0_dp
436  IF (ss*sscale > gcutoff) THEN
437  CALL xpbe_hole_t_c_lr_lda_calc_1(e_tmp, e_rho(ip), e_ndrho(ip), &
438  my_rho, &
439  my_ndrho, sscale, sx, r, order)
440  ELSE
441  CALL xpbe_hole_t_c_lr_lda_calc_2(e_tmp, e_rho(ip), e_ndrho(ip), &
442  my_rho, &
443  my_ndrho, sscale, sx, r, order)
444  END IF
445  e_0(ip) = e_0(ip) + 0.5_dp*e_tmp
446 
447  END IF
448  END DO
449 
450 !$OMP END DO
451 
452  END SUBROUTINE xpbe_hole_t_c_lr_lsd_calc
453 
454 ! **************************************************************************************************
455 !> \brief low level routine that calculates the energy derivatives in one point
456 !> \param e_0 derivatives of the energy on the grid
457 !> \param e_rho derivatives of the energy on the grid
458 !> \param e_ndrho derivatives of the energy on the grid
459 !> \param rho value of density on the grid
460 !> \param ndrho value of gradient on the grid
461 !> \param sscale functional parameters
462 !> \param sx functional parameters
463 !> \param R functional parameters
464 !> \param order order of the derivatives
465 !> \par History
466 !> 01.2009 created [Manuel Guidon]
467 !> \author Manuel Guidon
468 ! **************************************************************************************************
469  ELEMENTAL SUBROUTINE xpbe_hole_t_c_lr_lda_calc_1(e_0, e_rho, e_ndrho, &
470  rho, ndrho, sscale, sx, R, order)
471  REAL(kind=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho
472  REAL(kind=dp), INTENT(IN) :: rho, ndrho, sscale, sx, r
473  INTEGER, INTENT(IN) :: order
474 
475  REAL(kind=dp) :: p, pndrho, prho, q, qndrho, qrho, sndrho, srho, ssval, t1, t103, t104, &
476  t105, t106, t107, t108, t109, t11, t113, t115, t116, t117, t118, t119, t12, t121, t123, &
477  t125, t126, t129, t13, t133, t136, t14, t140, t142, t143, t144, t147, t150, t152, t155, &
478  t156, t159, t162, t163, t164, t167, t17, t170, t173, t175, t176, t177, t178, t181, t183, &
479  t187, t19, t190, t194, t195, t196, t199, t2, t202, t209, t21, t214, t216, t217, t218, &
480  t22, t222, t223, t224, t227, t228, t229, t232, t234, t235, t236, t237, t240, t244, t248, &
481  t25, t251, t255, t256, t257, t261, t262, t265, t269, t27, t273, t276
482  REAL(kind=dp) :: t279, t280, t281, t29, t292, t3, t308, t31, t312, t314, t318, t320, t322, &
483  t33, t331, t334, t339, t34, t341, t347, t349, t35, t359, t362, t368, t37, t370, t373, &
484  t38, t382, t385, t39, t393, t4, t40, t401, t402, t404, t405, t409, t41, t42, t420, t421, &
485  t424, t425, t426, t429, t430, t431, t437, t44, t445, t446, t45, t451, t46, t473, t475, &
486  t479, t48, t484, t497, t498, t5, t50, t503, t504, t51, t517, t52, t523, t524, t527, t53, &
487  t533, t534, t54, t541, t56, t568, t57, t58, t581, t590, t6, t60, t603, t604, t61, t611, &
488  t612, t64, t640, t65, t653, t667, t675, t677, t69, t7, t71, t716
489  REAL(kind=dp) :: t717, t720, t723, t73, t739, t758, t762, t77, t8, t800, t803, t808, t85, &
490  t86, t87, t88, t91, t92, t93, t94, t95, t98, t99
491 
492  IF (order >= 0) THEN
493  t1 = 3**(0.1e1_dp/0.3e1_dp)
494  t2 = t1**2
495  t3 = ndrho*t2
496  t4 = pi**2
497  t5 = t4*rho
498  t6 = t5**(0.1e1_dp/0.3e1_dp)
499  t7 = 0.1e1_dp/t6
500  t8 = 0.1e1_dp/rho
501  ssval = t3*t7*t8/0.6e1_dp
502  t11 = red(rho, ndrho)
503  t12 = t11**2
504  t13 = sscale**2
505  t14 = t12*t13
506  t17 = t12**2
507  t19 = t13**2
508  t21 = a1*t12*t13 + a2*t17*t19
509  t22 = t14*t21
510  t25 = t17*t11
511  t27 = t19*sscale
512  t29 = t17*t12
513  t31 = t19*t13
514  t33 = 1 + a3*t17*t19 + a4*t25*t27 + a5*t29*t31
515  t34 = 0.1e1_dp/t33
516  t35 = r**2
517  t37 = t6**2
518  t38 = t2*t37
519  t39 = t34*t35*t38
520  q = t22*t39
521  t40 = t21*t34
522  t41 = 0.1e1_dp/a
523  t42 = t40*t41
524  p = 0.9e1_dp/0.4e1_dp*t14*t42
525  t44 = rho**(0.1e1_dp/0.3e1_dp)
526  t45 = t44*rho
527  t46 = exei(p, q)
528  t48 = expint(1, q)
529  t50 = t14*t40
530  t51 = d + t50
531  t52 = t51*t35
532  t53 = t52*t38
533  t54 = expint(1, t53)
534  t56 = t51**2
535  t57 = t56*t51
536  t58 = 0.1e1_dp/t57
537  t60 = d*c
538  t61 = d**2
539  t64 = d*t21
540  t65 = t34*b
541  t69 = rootpi
542  t71 = f1*t21
543  t73 = t71*t34 + f2
544  t77 = c*(1 + t73*t12*t13)
545  t85 = t69*(15*e + 6*t77*t51 + 4*b*t56 + 8*a*t57)
546  t86 = sqrt(t51)
547  t87 = t86*t57
548  t88 = 0.1e1_dp/t87
549  t91 = sqrt(a)
550  t92 = pi*t91
551  t93 = exp(p)
552  t94 = t11*sscale
553  t95 = sqrt(t42)
554  t98 = erf(0.3e1_dp/0.2e1_dp*t94*t95)
555  t99 = 1 - t98
556  t103 = 0.3e1_dp/0.4e1_dp*pi + t85*t88/0.16e2_dp - 0.3e1_dp/0.4e1_dp*t92 &
557  *t93*t99
558  t104 = 0.1e1_dp/t69
559  t105 = t103*t104
560  t106 = 0.1e1_dp/t12
561  t107 = 0.1e1_dp/t13
562  t108 = t106*t107
563  t109 = t108*t87
564  t113 = t40*c + real(2*t64*t65, kind=dp) - 0.32e2_dp/0.15e2_dp*t105*t109 &
565  + t60*t73
566  t115 = t17*t19
567  t116 = t21**2
568  t117 = t33**2
569  t118 = 0.1e1_dp/t117
570  t119 = t116*t118
571  t121 = c*t73
572  t123 = t119*b + t40*t121
573  t125 = t35*t2
574  t126 = t61*c
575  t129 = e*t21
576  t133 = d*t103*t104
577  t136 = t34*c
578  t140 = real(2*t129*t34, kind=dp) - 0.32e2_dp/0.15e2_dp*t133*t109 + real(2 &
579  *t64*t136, kind=dp) + t126*t73
580  t142 = t105*t106
581  t143 = t107*t87
582  t144 = t143*t40
583  t147 = t136*t73
584  t150 = c*t116
585  t152 = -0.32e2_dp/0.15e2_dp*t142*t144 + real(2*t64*t147, kind=dp) + t150 &
586  *t118
587  t155 = t29*t31*c
588  t156 = t73*t116
589  t159 = t126 + 2*d*e + t14*t140 + t115*t152 + t155*t156* &
590  t118
591  t162 = t35**2
592  t163 = t162*t1
593  t164 = t6*t5
594  t167 = t61*t103*t104
595  t170 = t34*e
596  t173 = -0.16e2_dp/0.15e2_dp*t167*t109 + real(2*t64*t170, kind=dp)
597  t175 = t34*t103
598  t176 = t64*t175
599  t177 = t104*t106
600  t178 = t177*t143
601  t181 = e*t116
602  t183 = -0.32e2_dp/0.15e2_dp*t176*t178 + t181*t118
603  t187 = t104*t87*t119
604  t190 = t61*e + t14*t173 + t115*t183 - 0.16e2_dp/0.15e2_dp*t115 &
605  *t103*t187
606  t194 = 2*e + t60 + t61*b + t14*t113 + t115*t123 + t125*t37 &
607  *t159 + 3*t163*t164*t190
608  t195 = t58*t194
609  t196 = d*t35
610  t199 = exp(-t196*t38 - q)
611  t202 = -0.4e1_dp/0.9e1_dp*a*t46 + 0.4e1_dp/0.9e1_dp*a*t48 - 0.4e1_dp/ &
612  0.9e1_dp*a*t54 - 0.4e1_dp/0.9e1_dp*t195*t199
613  e_0 = e_0 + (t45*t202*clda)*sx
614  END IF
615  IF (order >= 1 .OR. order >= -1) THEN
616  t209 = rho**2
617  srho = -t3/t164*t8*t4/0.18e2_dp - t3*t7/t209/0.6e1_dp
618 
619  t214 = t2*t7
620  sndrho = t214*t8/0.6e1_dp
621  t216 = t11*t13
622  t217 = t216*t40
623  t218 = dsdrho(rho, ndrho)
624  t222 = 2*t217*t125*t37*t218
625  t223 = a1*t11
626  t224 = t13*t218
627  t227 = t12*t11
628  t228 = a2*t227
629  t229 = t19*t218
630  t232 = 2*t223*t224 + 4*t228*t229
631  t234 = t14*t232*t39
632  t235 = t21*t118
633  t236 = t14*t235
634  t237 = a3*t227
635  t240 = a4*t17
636  t244 = a5*t25
637  t248 = 4*t237*t229 + 5*t240*t27*t218 + 6*t244*t31*t218
638  t251 = t236*t125*t37*t248
639  t255 = 0.2e1_dp/0.3e1_dp*t50*t125*t7*t4
640  qrho = t222 + t234 - t251 + t255
641  t256 = t216*t21
642  t257 = t34*t41
643  t261 = t232*t34
644  t262 = t261*t41
645  t265 = t118*t41
646  prho = 0.9e1_dp/0.2e1_dp*t256*t257*t218 + 0.9e1_dp/0.4e1_dp*t14*t262 &
647  - 0.9e1_dp/0.4e1_dp*t22*t265*t248
648  t269 = dsdndrho(rho)
649  t273 = t13*t269
650  t276 = t19*t269
651  t279 = 2*t223*t273 + 4*t228*t276
652  t280 = t279*t34
653  t281 = t280*t41
654  t292 = 4*t237*t276 + 5*t240*t27*t269 + 6*t244*t31*t269
655  pndrho = 0.9e1_dp/0.2e1_dp*t256*t257*t269 + 0.9e1_dp/0.4e1_dp*t14* &
656  t281 - 0.9e1_dp/0.4e1_dp*t22*t265*t292
657  qndrho = 2*t217*t125*t37*t269 + t14*t279*t39 - t236*t125 &
658  *t37*t292
659  t308 = dexeirho(p, q, prho, qrho)
660  t312 = exp(-q)
661  t314 = t312*t106*t107
662  t318 = 0.1e1_dp/t35
663  t320 = 0.1e1_dp/t37
664  t322 = 0.1e1_dp/t21*t33*t318*t1*t320
665  t331 = 2*t216*t40*t218 + t14*t261 - t14*t235*t248
666  t334 = t214*t4
667  t339 = exp(-t53)
668  t341 = 0.1e1_dp/t51
669  t347 = t56**2
670  t349 = 0.1e1_dp/t347*t194
671  t359 = d*t232
672  t362 = t118*b
673  t368 = t118*t248
674  t370 = f1*t232*t34 - t71*t368
675  t373 = t73*t11
676  t382 = b*t51
677  t385 = a*t56
678  t393 = 0.1e1_dp/t86/t347
679  t401 = t92*t93
680  t402 = sqrt(0.31415926535897932385e1_dp)
681  t404 = exp(-p)
682  t405 = 0.1e1_dp/t402*t404
683  t409 = 0.1e1_dp/t95
684  t420 = real(t69*(6*c*(t370*t12*t13 + 2*t373*t224)*t51 &
685  + 6*t77*t331 + 8*t382*t331 + 24*t385*t331)*t88, kind=dp)/ &
686  0.16e2_dp - 0.7e1_dp/0.32e2_dp*real(t85, kind=dp)*real(t393, kind=dp)* &
687  REAL(t331, kind=dp) - 0.3e1_dp &
688  /0.4e1_dp*t92*prho*t93*t99 + 0.3e1_dp/0.2e1_dp*t401*t405 &
689  *(0.3e1_dp/0.2e1_dp*t218*sscale*t95 + 0.3e1_dp/0.4e1_dp*t94*t409 &
690  *(t262 - t235*t41*t248))
691  t421 = t420*t104
692  t424 = 0.1e1_dp/t227
693  t425 = t105*t424
694  t426 = t143*t218
695  t429 = t86*t56
696  t430 = t107*t429
697  t431 = t430*t331
698  t437 = t227*t19
699  t445 = 0.1e1_dp/t117/t33
700  t446 = t116*t445
701  t451 = t121*t248
702  t473 = t424*t107
703  t475 = t473*t87*t218
704  t479 = t108*t429*t331
705  t484 = t118*c
706  t497 = t105*t473
707  t498 = t87*t21
708  t503 = t105*t108
709  t504 = t429*t21
710  t517 = t64*t118
711  t523 = c*t21
712  t524 = t118*t232
713  t527 = t445*t248
714  t533 = t25*t31*c
715  t534 = t118*t218
716  t541 = t73*t21
717  t568 = t118*e
718  t581 = t64*t118*t103
719  t590 = t104*t424
720  t603 = t437*t105
721  t604 = t87*t116
722  t611 = t115*t105
723  t612 = t429*t116
724  e_rho = e_rho + (0.4e1_dp/0.3e1_dp*t44*t202*clda + t45*(-0.4e1_dp/0.9e1_dp* &
725  a*t308 - 0.4e1_dp/0.27e2_dp*a*qrho*t314*t322 + 0.4e1_dp/0.27e2_dp &
726  *a*(t331*t35*t38 + 0.2e1_dp/0.3e1_dp*t52*t334)*t339*t341 &
727  *t318*t1*t320 + 0.4e1_dp/0.3e1_dp*t349*t199*t331 - 0.4e1_dp &
728  /0.9e1_dp*t58*(real(2*t216*t113*t218, kind=dp) + t14*(t261*c &
729  - t235*c*real(t248, kind=dp) + real(2*t359*t65, kind=dp) - real(2*t64*t362 &
730  *t248, kind=dp) - 0.32e2_dp/0.15e2_dp*t421*t109 + 0.64e2_dp/0.15e2_dp*t425 &
731  *t426 - 0.112e3_dp/0.15e2_dp*t142*t431 + t60*t370) + real(4* &
732  t437*t123*t218, kind=dp) + t115*(0.2e1_dp*t235*b*t232 - 0.2e1_dp*t446 &
733  *b*real(t248, kind=dp) + t261*t121 - t235*t451 + t40*c*t370) &
734  + 0.2e1_dp/0.3e1_dp*t125*t7*t159*t4 + t125*t37*(real(2* &
735  t216*t140*t218, kind=dp) + t14*(0.2e1_dp*e*t232*t34 - real(2*t129 &
736  *t368, kind=dp) - 0.32e2_dp/0.15e2_dp*d*t420*t104*t109 + 0.64e2_dp/0.15e2_dp &
737  *t133*t475 - 0.112e3_dp/0.15e2_dp*t133*t479 + real(2*t359 &
738  *t136, kind=dp) - real(2*t64*t484*t248, kind=dp) + t126*t370) + real(4* &
739  t437*t152*t218, kind=dp) + t115*(-0.32e2_dp/0.15e2_dp*t421*t106*t144 &
740  + 0.64e2_dp/0.15e2_dp*t497*t498*t34*real(t218, kind=dp) - 0.112e3_dp/0.15e2_dp &
741  *t503*t504*t34*t331 - 0.32e2_dp/0.15e2_dp*t142*t143* &
742  t261 + 0.32e2_dp/0.15e2_dp*t503*t498*real(t368, kind=dp) + real(2*t359 &
743  *t147, kind=dp) - 0.2e1_dp*t517*t451 + 0.2e1_dp*real(t64, kind=dp)*real(t136, kind=dp)* &
744  t370 + real(2*t523*t524, kind=dp) - real(2*t150*t527, kind=dp)) + real(6*t533 &
745  *t156*t534, kind=dp) + t155*t370*t116*t118 + 0.2e1_dp*t155*t541 &
746  *real(t524, kind=dp) - 0.2e1_dp*t155*real(t156, kind=dp)*real(t527, kind=dp)) + 0.4e1_dp &
747  *t163*t6*t190*t4 + 0.3e1_dp*t163*t164*(real(2*t216*t173 &
748  *t218, kind=dp) + t14*(-0.16e2_dp/0.15e2_dp*t61*t420*t104*t109 + &
749  0.32e2_dp/0.15e2_dp*t167*t475 - 0.56e2_dp/0.15e2_dp*t167*t479 + real(2 &
750  *t359*t170, kind=dp) - real(2*t64*t568*t248, kind=dp)) + real(4*t437 &
751  *t183*t218, kind=dp) + t115*(-0.32e2_dp/0.15e2_dp*real(t359, kind=dp)*real(t175, kind=dp) &
752  *real(t178, kind=dp) + 0.32e2_dp/0.15e2_dp*t581*t177*t143*real(t248, kind=dp) &
753  - 0.32e2_dp/0.15e2_dp*real(t64, kind=dp)*t34*t420*real(t178, kind=dp) + 0.64e2_dp &
754  /0.15e2_dp*t176*t590*t426 - 0.112e3_dp/0.15e2_dp*t176*t177* &
755  t431 + real(2*t129*t524, kind=dp) - real(2*t181*t527, kind=dp)) - 0.64e2_dp/ &
756  0.15e2_dp*real(t603, kind=dp)*real(t604, kind=dp)*real(t534, kind=dp) - 0.16e2_dp/0.15e2_dp* &
757  t115*t420*t187 - 0.56e2_dp/0.15e2_dp*t611*t612*t118*t331 - &
758  0.32e2_dp/0.15e2_dp*t611*t498*real(t524, kind=dp) + 0.32e2_dp/0.15e2_dp*t611 &
759  *real(t604, kind=dp)*real(t527, kind=dp)))*t199 - 0.4e1_dp/0.9e1_dp*t195*(-0.2e1_dp &
760  /0.3e1_dp*t196*t334 - t222 - t234 + t251 - t255)*t199)* &
761  clda)*sx
762  t640 = dexeindrho(p, q, pndrho, qndrho)
763  t653 = 2*t216*t40*t269 + t14*t280 - t14*t235*t292
764  t667 = d*t279
765  t675 = t118*t292
766  t677 = f1*t279*t34 - t71*t675
767  t716 = real(t69*(6*c*(t677*t12*t13 + 2*t373*t273)*t51 &
768  + 6*t77*t653 + 8*t382*t653 + 24*t385*t653)*t88, kind=dp)/ &
769  0.16e2_dp - 0.7e1_dp/0.32e2_dp*real(t85, kind=dp)*real(t393, kind=dp)* &
770  REAL(t653, kind=dp) - 0.3e1_dp &
771  /0.4e1_dp*t92*pndrho*t93*t99 + 0.3e1_dp/0.2e1_dp*t401*t405 &
772  *(0.3e1_dp/0.2e1_dp*t269*sscale*t95 + 0.3e1_dp/0.4e1_dp*t94* &
773  t409*(t281 - t235*t41*t292))
774  t717 = t716*t104
775  t720 = t143*t269
776  t723 = t430*t653
777  t739 = t121*t292
778  t758 = t473*t87*t269
779  t762 = t108*t429*t653
780  t800 = t118*t279
781  t803 = t445*t292
782  t808 = t118*t269
783  e_ndrho = e_ndrho + (t45*(-0.4e1_dp/0.9e1_dp*a*t640 - 0.4e1_dp/0.27e2_dp*a*qndrho &
784  *t314*t322 + 0.4e1_dp/0.9e1_dp*a*t653*t339*t341 + 0.4e1_dp &
785  /0.3e1_dp*t349*t199*t653 - 0.4e1_dp/0.9e1_dp*t58*(real(2*t216 &
786  *t113*t269, kind=dp) + t14*(t280*c - t235*c*real(t292, kind=dp) + real(2 &
787  *t667*t65, kind=dp) - real(2*t64*t362*t292, kind=dp) - 0.32e2_dp/0.15e2_dp &
788  *t717*t109 + 0.64e2_dp/0.15e2_dp*t425*t720 - 0.112e3_dp/0.15e2_dp* &
789  t142*t723 + t60*t677) + real(4*t437*t123*t269, kind=dp) + t115* &
790  (0.2e1_dp*t235*b*t279 - 0.2e1_dp*t446*b*real(t292, kind=dp) + t280* &
791  t121 - t235*t739 + t40*c*t677) + t125*t37*(real(2*t216 &
792  *t140*t269, kind=dp) + t14*(0.2e1_dp*e*t279*t34 - real(2*t129* &
793  t675, kind=dp) - 0.32e2_dp/0.15e2_dp*d*t716*t104*t109 + 0.64e2_dp/0.15e2_dp &
794  *t133*t758 - 0.112e3_dp/0.15e2_dp*t133*t762 + real(2*t667* &
795  t136, kind=dp) - real(2*t64*t484*t292, kind=dp) + t126*t677) + real(4*t437 &
796  *t152*t269, kind=dp) + t115*(-0.32e2_dp/0.15e2_dp*t717*t106*t144 + &
797  0.64e2_dp/0.15e2_dp*t497*t498*t34*real(t269, kind=dp) - 0.112e3_dp/0.15e2_dp &
798  *t503*t504*t34*t653 - 0.32e2_dp/0.15e2_dp*t142*t143*t280 &
799  + 0.32e2_dp/0.15e2_dp*t503*t498*real(t675, kind=dp) + real(2*t667* &
800  t147, kind=dp) - 0.2e1_dp*t517*t739 + 0.2e1_dp*real(t64, kind=dp)*real(t136, kind=dp)*t677 &
801  + real(2*t523*t800, kind=dp) - real(2*t150*t803, kind=dp)) + real(6*t533 &
802  *t156*t808, kind=dp) + t155*t677*t116*t118 + 0.2e1_dp*t155*t541 &
803  *real(t800, kind=dp) - 0.2e1_dp*t155*real(t156, kind=dp)*real(t803, kind=dp)) + 0.3e1_dp*t163 &
804  *t164*(real(2*t216*t173*t269, kind=dp) + t14*(-0.16e2_dp/0.15e2_dp &
805  *t61*t716*t104*t109 + 0.32e2_dp/0.15e2_dp*t167*t758 - 0.56e2_dp &
806  /0.15e2_dp*t167*t762 + real(2*t667*t170, kind=dp) - real(2*t64 &
807  *t568*t292, kind=dp)) + real(4*t437*t183*t269, kind=dp) + t115*(-0.32e2_dp &
808  /0.15e2_dp*real(t667, kind=dp)*real(t175, kind=dp)*real(t178, kind=dp) + 0.32e2_dp/0.15e2_dp &
809  *t581*t177*t143*real(t292, kind=dp) - 0.32e2_dp/0.15e2_dp*real(t64, kind=dp)* &
810  t34*t716*real(t178, kind=dp) + 0.64e2_dp/0.15e2_dp*t176*t590*t720 - 0.112e3_dp &
811  /0.15e2_dp*t176*t177*t723 + real(2*t129*t800, kind=dp) - real(2 &
812  *t181*t803, kind=dp)) - 0.64e2_dp/0.15e2_dp*real(t603, kind=dp)*real(t604, kind=dp)* &
813  REAL(t808, kind=dp) - 0.16e2_dp/0.15e2_dp*t115*t716*t187 - 0.56e2_dp/0.15e2_dp &
814  *t611*t612*t118*t653 - 0.32e2_dp/0.15e2_dp*t611*t498*real(t800, kind=dp) &
815  + 0.32e2_dp/0.15e2_dp*t611*real(t604, kind=dp)*real(t803, kind=dp)))*t199 &
816  + 0.4e1_dp/0.9e1_dp*t195*qndrho*t199)*clda)*sx
817  END IF
818 
819  END SUBROUTINE xpbe_hole_t_c_lr_lda_calc_1
820 
821 ! **************************************************************************************************
822 !> \brief low level routine that calculates the energy derivatives in one point
823 !> \param e_0 derivatives of the energy on the grid
824 !> \param e_rho derivatives of the energy on the grid
825 !> \param e_ndrho derivatives of the energy on the grid
826 !> \param rho value of density on the grid
827 !> \param ndrho value of gradient on the grid
828 !> \param sscale functional parameters
829 !> \param sx functional parameters
830 !> \param R functional parameters
831 !> \param order order of the derivatives
832 !> \par History
833 !> 01.2009 created [Manuel Guidon]
834 !> \author Manuel Guidon
835 ! **************************************************************************************************
836  ELEMENTAL SUBROUTINE xpbe_hole_t_c_lr_lda_calc_2(e_0, e_rho, e_ndrho, &
837  rho, ndrho, sscale, sx, R, order)
838  REAL(kind=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho
839  REAL(kind=dp), INTENT(IN) :: rho, ndrho, sscale, sx, r
840  INTEGER, INTENT(IN) :: order
841 
842  REAL(kind=dp) :: p, pndrho, prho, q, qndrho, qrho, sndrho, srho, ssval, t1, t102, t106, t11, &
843  t110, t113, t115, t117, t118, t119, t12, t122, t125, t126, t127, t128, t13, t130, t133, &
844  t135, t138, t14, t140, t142, t143, t146, t150, t151, t152, t155, t158, t165, t17, t170, &
845  t172, t173, t174, t178, t179, t180, t183, t184, t185, t188, t19, t190, t191, t192, t193, &
846  t196, t197, t2, t200, t204, t207, t21, t211, t212, t213, t217, t22, t221, t225, t229, &
847  t232, t235, t236, t242, t248, t25, t264, t268, t27, t272, t276, t278, t280, t287, t289, &
848  t29, t292, t297, t299, t3, t305, t307, t31, t317, t320, t324
849  REAL(kind=dp) :: t327, t33, t330, t333, t334, t338, t34, t340, t344, t35, t352, t353, t358, &
850  t37, t38, t380, t39, t394, t398, t399, t4, t40, t401, t406, t407, t408, t41, t415, t435, &
851  t44, t45, t454, t46, t461, t48, t485, t496, t498, t5, t50, t51, t512, t52, t524, t525, &
852  t529, t53, t531, t54, t545, t56, t579, t58, t581, t586, t6, t60, t61, t64, t65, t7, t74, &
853  t75, t77, t79, t8, t81, t83, t84, t85, t86, t89, t91, t93, t94, t95, t97
854 
855  IF (order >= 0) THEN
856  t1 = 3**(0.1e1_dp/0.3e1_dp)
857  t2 = t1**2
858  t3 = ndrho*t2
859  t4 = pi**2
860  t5 = t4*rho
861  t6 = t5**(0.1e1_dp/0.3e1_dp)
862  t7 = 0.1e1_dp/t6
863  t8 = 0.1e1_dp/rho
864  ssval = t3*t7*t8/0.6e1_dp
865 
866  t11 = red(rho, ndrho)
867  t12 = t11**2
868  t13 = sscale**2
869  t14 = t12*t13
870  t17 = t12**2
871  t19 = t13**2
872  t21 = a1*t12*t13 + a2*t17*t19
873  t22 = t14*t21
874  t25 = t17*t11
875  t27 = t19*sscale
876  t29 = t17*t12
877  t31 = t19*t13
878  t33 = 1 + a3*t17*t19 + a4*t25*t27 + a5*t29*t31
879  t34 = 0.1e1_dp/t33
880  t35 = r**2
881  t37 = t6**2
882  t38 = t2*t37
883  t39 = t34*t35*t38
884  q = t22*t39
885  t40 = t21*t34
886  t41 = 0.1e1_dp/a
887  p = 0.9e1_dp/0.4e1_dp*t14*t40*t41
888  t44 = rho**(0.1e1_dp/0.3e1_dp)
889  t45 = t44*rho
890  t46 = exei(p, q)
891  t48 = expint(1, q)
892  t50 = t14*t40
893  t51 = d + t50
894  t52 = t51*t35
895  t53 = t52*t38
896  t54 = expint(1, t53)
897  t56 = t51**2
898  t58 = 0.1e1_dp/t56/t51
899  t60 = d*c
900  t61 = d**2
901  t64 = d*t21
902  t65 = t34*b
903  t74 = g1 + g2*t12*t13 + g3*t17*t19 + g4*t25*t27
904  t75 = e*t74
905  t77 = f1*t21
906  t79 = t77*t34 + f2
907  t81 = t40*c + 2*t64*t65 + 2*t75 + t60*t79
908  t83 = t17*t19
909  t84 = t21**2
910  t85 = t33**2
911  t86 = 0.1e1_dp/t85
912  t89 = c*t79
913  t91 = t84*t86*b + t40*t89
914  t93 = t35*t2
915  t94 = t61*c
916  t95 = d*e
917  t97 = e*t21
918  t102 = t34*c
919  t106 = 2*t97*t34 + 2*t95*t74 + 2*t64*t102 + t94*t79
920  t110 = t102*t79
921  t113 = c*t84
922  t115 = 2*t75*t40 + 2*t64*t110 + t113*t86
923  t117 = t29*t31
924  t118 = t117*c
925  t119 = t79*t84
926  t122 = t94 + 2*t95 + t14*t106 + t83*t115 + t118*t119*t86
927  t125 = t35**2
928  t126 = t125*t1
929  t127 = t6*t5
930  t128 = t61*e
931  t130 = t34*e
932  t133 = t128*t74 + 2*t64*t130
933  t135 = t130*t74
934  t138 = e*t84
935  t140 = 2*t64*t135 + t138*t86
936  t142 = t117*e
937  t143 = t74*t84
938  t146 = t128 + t14*t133 + t83*t140 + t142*t143*t86
939  t150 = 2*e + t60 + t61*b + t14*t81 + t83*t91 + t93*t37* &
940  t122 + 3*t126*t127*t146
941  t151 = t58*t150
942  t152 = d*t35
943  t155 = exp(-t152*t38 - q)
944  t158 = -0.4e1_dp/0.9e1_dp*a*t46 + 0.4e1_dp/0.9e1_dp*a*t48 - 0.4e1_dp/ &
945  0.9e1_dp*a*t54 - 0.4e1_dp/0.9e1_dp*t151*t155
946  e_0 = e_0 + (t45*t158*clda)*sx
947  END IF
948  IF (order >= 1 .OR. order == -1) THEN
949  t165 = rho**2
950  srho = -t3/t127*t8*t4/0.18e2_dp - t3*t7/t165/0.6e1_dp
951  t170 = t2*t7
952  sndrho = t170*t8/0.6e1_dp
953  t172 = t11*t13
954  t173 = t172*t40
955  t174 = dsdrho(rho, ndrho)
956  t178 = 2*t173*t93*t37*t174
957  t179 = a1*t11
958  t180 = t13*t174
959  t183 = t12*t11
960  t184 = a2*t183
961  t185 = t19*t174
962  t188 = 2*t179*t180 + 4*t184*t185
963  t190 = t14*t188*t39
964  t191 = t21*t86
965  t192 = t14*t191
966  t193 = a3*t183
967  t196 = a4*t17
968  t197 = t27*t174
969  t200 = a5*t25
970  t204 = 4*t193*t185 + 5*t196*t197 + 6*t200*t31*t174
971  t207 = t192*t93*t37*t204
972  t211 = 0.2e1_dp/0.3e1_dp*t50*t93*t7*t4
973  qrho = t178 + t190 - t207 + t211
974  t212 = t172*t21
975  t213 = t34*t41
976  t217 = t188*t34
977  t221 = t86*t41
978  prho = 0.9e1_dp/0.2e1_dp*t212*t213*t174 + 0.9e1_dp/0.4e1_dp*t14*t217 &
979  *t41 - 0.9e1_dp/0.4e1_dp*t22*t221*t204
980  t225 = dsdndrho(rho)
981  t229 = t13*t225
982  t232 = t19*t225
983  t235 = 2*t179*t229 + 4*t184*t232
984  t236 = t235*t34
985  t242 = t27*t225
986  t248 = 4*t193*t232 + 5*t196*t242 + 6*t200*t31*t225
987  pndrho = 0.9e1_dp/0.2e1_dp*t212*t213*t225 + 0.9e1_dp/0.4e1_dp*t14* &
988  t236*t41 - 0.9e1_dp/0.4e1_dp*t22*t221*t248
989  qndrho = 2*t173*t93*t37*t225 + t14*t235*t39 - t192*t93 &
990  *t37*t248
991  t264 = dexeirho(p, q, prho, qrho)
992  t268 = exp(-q)
993  t272 = t268/t12/t13
994  t276 = 0.1e1_dp/t35
995  t278 = 0.1e1_dp/t37
996  t280 = 0.1e1_dp/t21*t33*t276*t1*t278
997  t287 = t191*t204
998  t289 = 2*t172*t40*t174 + t14*t217 - t14*t287
999  t292 = t170*t4
1000  t297 = exp(-t53)
1001  t299 = 0.1e1_dp/t51
1002  t305 = t56**2
1003  t307 = 0.1e1_dp/t305*t150
1004  t317 = d*t188
1005  t320 = t86*b
1006  t324 = g2*t11
1007  t327 = g3*t183
1008  t330 = g4*t17
1009  t333 = 2*t324*t180 + 4*t327*t185 + 5*t330*t197
1010  t334 = e*t333
1011  t338 = t86*t204
1012  t340 = f1*t188*t34 - t77*t338
1013  t344 = t183*t19
1014  t352 = 0.1e1_dp/t85/t33
1015  t353 = t84*t352
1016  t358 = t89*t204
1017  t380 = t86*c
1018  t394 = t64*t86
1019  t398 = c*t21
1020  t399 = t86*t188
1021  t401 = t352*t204
1022  t406 = t25*t31
1023  t407 = t406*c
1024  t408 = t86*t174
1025  t415 = t79*t21
1026  t435 = t86*e
1027  t454 = t406*e
1028  t461 = t74*t21
1029  e_rho = e_rho + (0.4e1_dp/0.3e1_dp*t44*t158*clda + t45*(-0.4e1_dp/0.9e1_dp* &
1030  a*t264 - 0.4e1_dp/0.27e2_dp*a*qrho*t272*t280 + 0.4e1_dp/0.27e2_dp &
1031  *a*(t289*t35*t38 + 0.2e1_dp/0.3e1_dp*t52*t292)*t297*t299 &
1032  *t276*t1*t278 + 0.4e1_dp/0.3e1_dp*t307*t155*t289 - 0.4e1_dp &
1033  /0.9e1_dp*t58*(real(2*t172*t81*t174, kind=dp) + real(t14*(t217 &
1034  *c - t191*c*t204 + 2*t317*t65 - 2*t64*t320*t204 + 2 &
1035  *t334 + t60*t340), kind=dp) + real(4*t344*t91*t174, kind=dp) + real(t83* &
1036  (2*t191*b*t188 - 2*t353*b*t204 + t217*t89 - t191*t358 &
1037  + t40*c*t340), kind=dp) + 0.2e1_dp/0.3e1_dp*t93*t7*t122*t4 + t93 &
1038  *t37*real(2*t172*t106*t174 + t14*(2*e*t188*t34 &
1039  - 2*t97*t338 + 2*t95*t333 + 2*t317*t102 - 2*t64*t380 &
1040  *t204 + t94*t340) + 4*t344*t115*t174 + t83*(2*t334 &
1041  *t40 + 2*t75*t217 - 2*t75*t287 + 2*t317*t110 - 2*t394 &
1042  *t358 + 2*t64*t102*t340 + 2*t398*t399 - 2*t113* &
1043  t401) + 6*t407*t119*t408 + t118*t340*t84*t86 + 2*t118 &
1044  *t415*t399 - 2*t118*t119*t401, kind=dp) + 0.4e1_dp*t126*t6*t146 &
1045  *t4 + 0.3e1_dp*t126*t127*real(2*t172*t133*t174 + t14 &
1046  *(t128*t333 + 2*t317*t130 - 2*t64*t435*t204) + 4*t344 &
1047  *t140*t174 + t83*(2*t317*t135 - 2*t394*t75*t204 &
1048  + 2*t64*t130*t333 + 2*t97*t399 - 2*t138*t401) + 6* &
1049  t454*t143*t408 + t142*t333*t84*t86 + 2*t142*t461*t399 &
1050  - 2*t142*t143*t401, kind=dp))*t155 - 0.4e1_dp/0.9e1_dp*t151*(-0.2e1_dp &
1051  /0.3e1_dp*t152*t292 - t178 - t190 + t207 - t211)*t155)* &
1052  clda)*sx
1053  t485 = dexeindrho(p, q, pndrho, qndrho)
1054  t496 = t191*t248
1055  t498 = 2*t172*t40*t225 + t14*t236 - t14*t496
1056  t512 = d*t235
1057  t524 = 2*t324*t229 + 4*t327*t232 + 5*t330*t242
1058  t525 = e*t524
1059  t529 = t86*t248
1060  t531 = f1*t235*t34 - t77*t529
1061  t545 = t89*t248
1062  t579 = t86*t235
1063  t581 = t352*t248
1064  t586 = t86*t225
1065  e_ndrho = e_ndrho + (t45*(-0.4e1_dp/0.9e1_dp*a*t485 - 0.4e1_dp/0.27e2_dp*a*qndrho &
1066  *t272*t280 + 0.4e1_dp/0.9e1_dp*a*t498*t297*t299 + 0.4e1_dp &
1067  /0.3e1_dp*t307*t155*t498 - 0.4e1_dp/0.9e1_dp*t58*real(2*t172 &
1068  *t81*t225 + t14*(t236*c - t191*c*t248 + 2*t512*t65 &
1069  - 2*t64*t320*t248 + 2*t525 + t60*t531) + 4*t344*t91 &
1070  *t225 + t83*(2*t191*b*t235 - 2*t353*b*t248 + t236 &
1071  *t89 - t191*t545 + t40*c*t531) + t93*t37*(2*t172* &
1072  t106*t225 + t14*(2*e*t235*t34 - 2*t97*t529 + 2*t95 &
1073  *t524 + 2*t512*t102 - 2*t64*t380*t248 + t94*t531) + &
1074  4*t344*t115*t225 + t83*(2*t525*t40 + 2*t75*t236 - &
1075  2*t75*t496 + 2*t512*t110 - 2*t394*t545 + 2*t64*t102 &
1076  *t531 + 2*t398*t579 - 2*t113*t581) + 6*t407*t119* &
1077  t586 + t118*t531*t84*t86 + 2*t118*t415*t579 - 2*t118 &
1078  *t119*t581) + 3*t126*t127*(2*t172*t133*t225 + t14 &
1079  *(t128*t524 + 2*t512*t130 - 2*t64*t435*t248) + 4*t344 &
1080  *t140*t225 + t83*(2*t512*t135 - 2*t394*t75*t248 &
1081  + 2*t64*t130*t524 + 2*t97*t579 - 2*t138*t581) + 6* &
1082  t454*t143*t586 + t142*t524*t84*t86 + 2*t142*t461*t579 &
1083  - 2*t142*t143*t581), kind=dp)*t155 + 0.4e1_dp/0.9e1_dp*t151*qndrho &
1084  *t155)*clda)*sx
1085  END IF
1086 
1087  END SUBROUTINE xpbe_hole_t_c_lr_lda_calc_2
1088 
1089 ! **************************************************************************************************
1090 !> \brief These functions evaluate products exp(x)*Ei(x) and pi*exp(x)*erfc(sqrt(x)),
1091 !> as well as their derivatives with respect to various combinations of
1092 !> rho and norm_drho.
1093 !> \param P = 9/4*s^2*H/A
1094 !> \param Q = s^2*H*R^2*kf^2
1095 !> \return ...
1096 !> \par History
1097 !> 01.2009 created [Manuel Guidon]
1098 !> \author Manuel Guidon
1099 !> \note
1100 !> - In order to avoid numerical instabilities, these routines use Taylor-
1101 !> expansions for the above core-products for large arguments.
1102 !> - red calculates the reduced gradient in a save way (i.e. using a fixed
1103 !> cutoff EPS_RHO)
1104 ! **************************************************************************************************
1105  ELEMENTAL FUNCTION exei(P, Q)
1106  REAL(dp), INTENT(IN) :: p, q
1107  REAL(dp) :: exei
1108 
1109  REAL(dp) :: q2, q3, q4, tmp
1110 
1111  exei = 0.0_dp
1112  IF (p < expcutoff) THEN
1113  !Use exact product
1114  IF (p + q < 0.5_dp) THEN
1115  tmp = -euler - log(p + q) + p + q
1116  tmp = tmp - 0.25_dp*(p + q)**2 + 1.0_dp/18.0_dp*(p + q)**3 - 1.0_dp/96.0_dp*(p + q)**4
1117  tmp = tmp + 1.0_dp/600.0_dp*(p + q)**5
1118  exei = exp(p)*tmp
1119  ELSE
1120  exei = exp(p)*expint(1, q + p)
1121  END IF
1122  ELSE
1123  !Use approximation
1124  tmp = 1.0_dp/p
1125  ! *** 1st order
1126  exei = tmp
1127  ! *** 2nd order
1128  tmp = tmp/p
1129  exei = exei - (q + 1.0_dp)*tmp
1130  ! *** 3rd order
1131  tmp = tmp/p
1132  q2 = q*q
1133  exei = exei + (2.0_dp*q + q2 + 2.0_dp)*tmp
1134  ! *** 4th order
1135  tmp = tmp/p
1136  q3 = q2*q
1137  exei = exei - (3.0_dp*q2 + 6.0_dp*q + q3 + 6.0_dp)*tmp
1138  ! *** 5th order
1139  tmp = tmp/p
1140  q4 = q3*q
1141  exei = exei + (24.0_dp + 4.0_dp*q3 + q4 + 12.0_dp*q2 + 24.0_dp*q)*tmp
1142 
1143  ! *** scaling factor
1144  exei = exp(-q)*exei
1145  END IF
1146  END FUNCTION exei
1147 
1148 ! **************************************************************************************************
1149 !> \brief ...
1150 !> \param P ...
1151 !> \param Q ...
1152 !> \param dPrho ...
1153 !> \param dQrho ...
1154 !> \return ...
1155 ! **************************************************************************************************
1156  ELEMENTAL FUNCTION dexeirho(P, Q, dPrho, dQrho)
1157  REAL(dp), INTENT(IN) :: p, q, dprho, dqrho
1158  REAL(dp) :: dexeirho
1159 
1160  dexeirho = dprho*exei(p, q) - (dprho + dqrho)/(p + q)*exp(-q)
1161  END FUNCTION dexeirho
1162 
1163 ! **************************************************************************************************
1164 !> \brief ...
1165 !> \param P ...
1166 !> \param Q ...
1167 !> \param dPndrho ...
1168 !> \param dQndrho ...
1169 !> \return ...
1170 ! **************************************************************************************************
1171  ELEMENTAL FUNCTION dexeindrho(P, Q, dPndrho, dQndrho)
1172  REAL(dp), INTENT(IN) :: p, q, dpndrho, dqndrho
1173  REAL(dp) :: dexeindrho
1174 
1175  dexeindrho = dpndrho*exei(p, q) - (dpndrho + dqndrho)/(p + q)*exp(-q)
1176  END FUNCTION dexeindrho
1177 
1178 ! **************************************************************************************************
1179 !> \brief ...
1180 !> \param rho ...
1181 !> \param ndrho ...
1182 !> \return ...
1183 ! **************************************************************************************************
1184  ELEMENTAL FUNCTION red(rho, ndrho)
1185  REAL(dp), INTENT(IN) :: rho, ndrho
1186  REAL(dp) :: red
1187 
1188  red = 1.0_dp/6.0_dp*ndrho*3.0_dp**(2.0_dp/3.0_dp)
1189  red = red/(pi**(2.0_dp/3.0_dp))
1190  red = red*max(1.0_dp/rho**(4.0_dp/3.0_dp), eps_rho)
1191  END FUNCTION red
1192 
1193 ! **************************************************************************************************
1194 !> \brief ...
1195 !> \param rho ...
1196 !> \param ndrho ...
1197 !> \return ...
1198 ! **************************************************************************************************
1199  ELEMENTAL FUNCTION dsdrho(rho, ndrho)
1200  REAL(dp), INTENT(IN) :: rho, ndrho
1201  REAL(dp) :: dsdrho
1202 
1203  dsdrho = -2.0_dp/9.0_dp*ndrho*3.0**(2.0_dp/3.0_dp)
1204  dsdrho = dsdrho/(pi**(2.0_dp/3.0_dp))
1205  dsdrho = dsdrho*max(1.0_dp/rho**(7.0_dp/3.0_dp), eps_rho)
1206  END FUNCTION dsdrho
1207 
1208 ! **************************************************************************************************
1209 !> \brief ...
1210 !> \param rho ...
1211 !> \return ...
1212 ! **************************************************************************************************
1213  ELEMENTAL FUNCTION dsdndrho(rho)
1214  REAL(dp), INTENT(IN) :: rho
1215  REAL(dp) :: dsdndrho
1216 
1217  dsdndrho = 1.0_dp/6.0_dp*3.0_dp**(2.0_dp/3.0_dp)
1218  dsdndrho = dsdndrho/(pi**(2.0_dp/3.0_dp))
1219  dsdndrho = dsdndrho*max(1.0_dp/rho**(4.0_dp/3.0_dp), eps_rho)
1220  END FUNCTION dsdndrho
1221 
1222 ! **************************************************************************************************
1223 !> \brief computes the exponential integral
1224 !> En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,1,..
1225 !> Note: Ei(-x) = -E1(x)
1226 !> \param n ...
1227 !> \param x ...
1228 !> \return ...
1229 !> \par History
1230 !> 05.2007 Created
1231 !> \author Manuel Guidon (adapted from Numerical recipies, cleaner version of mathlib)
1232 ! **************************************************************************************************
1233  ELEMENTAL FUNCTION expint(n, x)
1234  INTEGER, INTENT(IN) :: n
1235  REAL(dp), INTENT(IN) :: x
1236  REAL(dp) :: expint
1237 
1238  INTEGER, PARAMETER :: maxit = 100
1239  REAL(dp), PARAMETER :: eps = 6.e-14_dp, euler = 0.5772156649015328606065120_dp, &
1240  fpmin = tiny(0.0_dp)
1241 
1242  INTEGER :: i, ii, nm1
1243  REAL(dp) :: a, b, c, d, del, fact, h, psi
1244 
1245  nm1 = n - 1
1246 
1247  IF (n .LT. 0 .OR. x .LT. 0.0_dp .OR. (x .EQ. 0.0_dp .AND. (n .EQ. 0 .OR. n .EQ. 1))) THEN
1248  expint = huge(1.0_dp)
1249  ELSE IF (n .EQ. 0) THEN !Special case.
1250  expint = exp(-x)/x
1251  ELSE IF (x .EQ. 0.0_dp) THEN !Another special case.
1252  expint = 1.0_dp/nm1
1253  ELSE IF (x .GT. 1.0_dp) THEN !Lentz's algorithm (5.2).
1254  b = x + n
1255  c = 1.0_dp/fpmin
1256  d = 1.0_dp/b
1257  h = d
1258  DO i = 1, maxit
1259  a = -i*(nm1 + i)
1260  b = b + 2.0_dp
1261  d = 1.0_dp/(a*d + b)
1262  c = b + a/c
1263  del = c*d
1264  h = h*del
1265  IF (abs(del - 1.0_dp) .LT. eps .OR. i == maxit) THEN
1266  expint = h*exp(-x)
1267  RETURN
1268  END IF
1269  END DO
1270  ELSE !Evaluate series.
1271  IF (nm1 .NE. 0) THEN !Set first term.
1272  expint = 1.0_dp/nm1
1273  ELSE
1274  expint = -log(x) - euler
1275  END IF
1276  fact = 1.0_dp
1277  DO i = 1, maxit
1278  fact = -fact*x/i
1279  IF (i .NE. nm1) THEN
1280  del = -fact/(i - nm1)
1281  ELSE
1282  psi = -euler !Compute I(n).
1283  DO ii = 1, nm1
1284  psi = psi + 1.0_dp/ii
1285  END DO
1286  del = fact*(-log(x) + psi)
1287  END IF
1288  expint = expint + del
1289  IF (abs(del) .LT. abs(expint)*eps) RETURN
1290  END DO
1291  END IF
1292  RETURN
1293  END FUNCTION expint
1294 
1295 END MODULE xc_xpbe_hole_t_c_lr
1296 
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
mathconstants::euler
real(kind=dp), parameter, public euler
Definition: mathconstants.F:119
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition: mathconstants.F:114
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_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_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_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_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_info
subroutine, public xpbe_hole_t_c_lr_lda_info(reference, shortform, needs, max_deriv)
returns various information on the functional
Definition: xc_xpbe_hole_t_c_lr.F:91
xc_xpbe_hole_t_c_lr::xpbe_hole_t_c_lr_lda_eval
subroutine, public xpbe_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
evaluates the pbe-hole exchange in a truncated coulomb potential
Definition: xc_xpbe_hole_t_c_lr.F:154
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_lsd_eval
subroutine, public xpbe_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params)
evaluates the pbe-hole exchange in a truncated coulomb potential
Definition: xc_xpbe_hole_t_c_lr.F:304
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_lsd_info
subroutine, public xpbe_hole_t_c_lr_lsd_info(reference, shortform, needs, max_deriv)
returns various information on the functional
Definition: xc_xpbe_hole_t_c_lr.F:121