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