Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 1626 : 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 1626 : IF (PRESENT(reference)) THEN
96 2 : reference = "{LDA version}"
97 : END IF
98 1626 : IF (PRESENT(shortform)) THEN
99 2 : shortform = "{LDA}"
100 : END IF
101 1626 : IF (PRESENT(needs)) THEN
102 1624 : needs%rho = .TRUE.
103 1624 : needs%norm_drho = .TRUE.
104 : END IF
105 1626 : IF (PRESENT(max_deriv)) max_deriv = 1
106 :
107 1626 : 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 792 : 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 792 : IF (PRESENT(reference)) THEN
126 2 : reference = "{LSD version}"
127 : END IF
128 792 : IF (PRESENT(shortform)) THEN
129 2 : shortform = "{LSD}"
130 : END IF
131 792 : IF (PRESENT(needs)) THEN
132 790 : needs%rho_spin = .TRUE.
133 790 : needs%norm_drho_spin = .TRUE.
134 : END IF
135 792 : IF (PRESENT(max_deriv)) max_deriv = 1
136 :
137 792 : 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 7008 : 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 1752 : POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
167 1752 : rho
168 : TYPE(xc_derivative_type), POINTER :: deriv
169 :
170 1752 : CALL timeset(routineN, handle)
171 :
172 1752 : CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
173 1752 : CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
174 :
175 : CALL xc_rho_set_get(rho_set, rho=rho, &
176 1752 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
177 1752 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
178 :
179 1752 : dummy => rho
180 :
181 1752 : e_0 => dummy
182 1752 : e_rho => dummy
183 1752 : e_ndrho => dummy
184 :
185 1752 : IF (order >= 0) THEN
186 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
187 1752 : allocate_deriv=.TRUE.)
188 1752 : CALL xc_derivative_get(deriv, deriv_data=e_0)
189 : END IF
190 1752 : IF (order >= 1 .OR. order == -1) THEN
191 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
192 1454 : allocate_deriv=.TRUE.)
193 1454 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
194 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
195 1454 : allocate_deriv=.TRUE.)
196 1454 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
197 : END IF
198 1752 : IF (order > 1 .OR. order < -1) THEN
199 0 : CPABORT("derivatives bigger than 1 not implemented")
200 : END IF
201 :
202 1752 : IF (R == 0.0_dp) THEN
203 0 : 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 1752 : !$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 1752 : CALL timestop(handle)
219 :
220 1752 : 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 1752 : 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 44020446 : my_rho = MAX(rho(ip), 0.0_dp)
257 44020446 : IF (my_rho > epsilon_rho) THEN
258 38464621 : 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 38464621 : sscale = 1.0_dp
261 38464621 : t1 = pi**2
262 38464621 : t2 = t1*my_rho
263 38464621 : t3 = t2**(0.1e1_dp/0.3e1_dp)
264 38464621 : t4 = 0.1e1_dp/t3
265 38464621 : t6 = my_ndrho*t4
266 38464621 : t7 = 0.1e1_dp/my_rho
267 38464621 : t8 = t7*sscale
268 38464621 : ss = 0.3466806371753173524216762e0_dp*t6*t8
269 38464621 : IF (ss > scutoff) THEN
270 24758866 : ss2 = ss*ss
271 24758866 : sscale = (smax*ss2 - sconst)/(ss2*ss)
272 : END IF
273 38464621 : 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 38439648 : 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 24973 : my_ndrho, sscale, sx, R, order)
281 : END IF
282 : END IF
283 : END DO
284 :
285 : !$OMP END DO
286 :
287 1752 : 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 3256 : 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 814 : POINTER :: dummy, e_0, e_ndrhoa, e_ndrhob, e_rhoa, &
317 814 : e_rhob, norm_drhoa, norm_drhob, rhoa, &
318 814 : rhob
319 : TYPE(xc_derivative_type), POINTER :: deriv
320 :
321 814 : CALL timeset(routineN, handle)
322 :
323 814 : CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
324 814 : 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 814 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, local_bounds=bo, rho_cutoff=epsilon_rho)
328 814 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
329 :
330 814 : dummy => rhoa
331 :
332 814 : e_0 => dummy
333 814 : e_rhoa => dummy
334 814 : e_rhob => dummy
335 814 : e_ndrhoa => dummy
336 814 : e_ndrhob => dummy
337 :
338 814 : IF (order >= 0) THEN
339 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
340 814 : allocate_deriv=.TRUE.)
341 814 : CALL xc_derivative_get(deriv, deriv_data=e_0)
342 : END IF
343 814 : IF (order >= 1 .OR. order == -1) THEN
344 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
345 628 : allocate_deriv=.TRUE.)
346 628 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
347 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
348 628 : allocate_deriv=.TRUE.)
349 628 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
350 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
351 628 : allocate_deriv=.TRUE.)
352 628 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
353 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
354 628 : allocate_deriv=.TRUE.)
355 628 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
356 : END IF
357 814 : IF (order > 1 .OR. order < -1) THEN
358 0 : CPABORT("derivatives bigger than 1 not implemented")
359 : END IF
360 :
361 814 : IF (R == 0.0_dp) THEN
362 0 : 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 814 : !$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 814 : CALL timestop(handle)
380 :
381 814 : 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 1628 : 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 27647946 : my_rho = MAX(2.0_dp*rho(ip), 0.0_dp)
418 27647946 : IF (my_rho > epsilon_rho) THEN
419 27540550 : 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 27540550 : sscale = 1.0_dp
423 27540550 : t1 = pi**2
424 27540550 : t2 = t1*my_rho
425 27540550 : t3 = t2**(0.1e1_dp/0.3e1_dp)
426 27540550 : t4 = 0.1e1_dp/t3
427 27540550 : t6 = my_ndrho*t4
428 27540550 : t7 = 0.1e1_dp/my_rho
429 27540550 : t8 = t7*sscale
430 27540550 : ss = 0.3466806371753173524216762e0_dp*t6*t8
431 27540550 : IF (ss > scutoff) THEN
432 12159382 : ss2 = ss*ss
433 12159382 : sscale = (smax*ss2 - sconst)/(ss2*ss)
434 : END IF
435 27540550 : e_tmp = 0.0_dp
436 27540550 : 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 27533842 : 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 6708 : my_ndrho, sscale, sx, R, order)
444 : END IF
445 27540550 : 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 1628 : 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 77657996 : 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 77657996 : IF (order >= 0) THEN
493 77657996 : t1 = 3**(0.1e1_dp/0.3e1_dp)
494 77657996 : t2 = t1**2
495 77657996 : t3 = ndrho*t2
496 77657996 : t4 = pi**2
497 77657996 : t5 = t4*rho
498 77657996 : t6 = t5**(0.1e1_dp/0.3e1_dp)
499 77657996 : t7 = 0.1e1_dp/t6
500 77657996 : t8 = 0.1e1_dp/rho
501 77657996 : ssval = t3*t7*t8/0.6e1_dp
502 77657996 : t11 = red(rho, ndrho)
503 77657996 : t12 = t11**2
504 77657996 : t13 = sscale**2
505 77657996 : t14 = t12*t13
506 77657996 : t17 = t12**2
507 77657996 : t19 = t13**2
508 77657996 : t21 = a1*t12*t13 + a2*t17*t19
509 77657996 : t22 = t14*t21
510 77657996 : t25 = t17*t11
511 77657996 : t27 = t19*sscale
512 77657996 : t29 = t17*t12
513 77657996 : t31 = t19*t13
514 77657996 : t33 = 1 + a3*t17*t19 + a4*t25*t27 + a5*t29*t31
515 77657996 : t34 = 0.1e1_dp/t33
516 77657996 : t35 = R**2
517 77657996 : t37 = t6**2
518 77657996 : t38 = t2*t37
519 77657996 : t39 = t34*t35*t38
520 77657996 : q = t22*t39
521 77657996 : t40 = t21*t34
522 77657996 : t41 = 0.1e1_dp/A
523 77657996 : t42 = t40*t41
524 77657996 : p = 0.9e1_dp/0.4e1_dp*t14*t42
525 77657996 : t44 = rho**(0.1e1_dp/0.3e1_dp)
526 77657996 : t45 = t44*rho
527 77657996 : t46 = exei(P, Q)
528 77657996 : t48 = expint(1, q)
529 77657996 : t50 = t14*t40
530 77657996 : t51 = D + t50
531 77657996 : t52 = t51*t35
532 77657996 : t53 = t52*t38
533 77657996 : t54 = expint(1, t53)
534 77657996 : t56 = t51**2
535 77657996 : t57 = t56*t51
536 77657996 : t58 = 0.1e1_dp/t57
537 77657996 : t60 = D*C
538 77657996 : t61 = D**2
539 77657996 : t64 = D*t21
540 77657996 : t65 = t34*B
541 77657996 : t69 = rootpi
542 77657996 : t71 = F1*t21
543 77657996 : t73 = t71*t34 + F2
544 77657996 : t77 = C*(1 + t73*t12*t13)
545 77657996 : t85 = t69*(15*E + 6*t77*t51 + 4*B*t56 + 8*A*t57)
546 77657996 : t86 = SQRT(t51)
547 77657996 : t87 = t86*t57
548 77657996 : t88 = 0.1e1_dp/t87
549 77657996 : t91 = SQRT(A)
550 77657996 : t92 = pi*t91
551 77657996 : t93 = EXP(p)
552 77657996 : t94 = t11*sscale
553 77657996 : t95 = SQRT(t42)
554 77657996 : t98 = erf(0.3e1_dp/0.2e1_dp*t94*t95)
555 77657996 : 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 77657996 : *t93*t99
558 77657996 : t104 = 0.1e1_dp/t69
559 77657996 : t105 = t103*t104
560 77657996 : t106 = 0.1e1_dp/t12
561 77657996 : t107 = 0.1e1_dp/t13
562 77657996 : t108 = t106*t107
563 77657996 : t109 = t108*t87
564 : t113 = t40*C + REAL(2*t64*t65, KIND=dp) - 0.32e2_dp/0.15e2_dp*t105*t109 &
565 77657996 : + t60*t73
566 77657996 : t115 = t17*t19
567 77657996 : t116 = t21**2
568 77657996 : t117 = t33**2
569 77657996 : t118 = 0.1e1_dp/t117
570 77657996 : t119 = t116*t118
571 77657996 : t121 = C*t73
572 77657996 : t123 = t119*B + t40*t121
573 77657996 : t125 = t35*t2
574 77657996 : t126 = t61*C
575 77657996 : t129 = E*t21
576 77657996 : t133 = D*t103*t104
577 77657996 : t136 = t34*C
578 : t140 = REAL(2*t129*t34, KIND=dp) - 0.32e2_dp/0.15e2_dp*t133*t109 + REAL(2 &
579 77657996 : *t64*t136, KIND=dp) + t126*t73
580 77657996 : t142 = t105*t106
581 77657996 : t143 = t107*t87
582 77657996 : t144 = t143*t40
583 77657996 : t147 = t136*t73
584 77657996 : t150 = C*t116
585 : t152 = -0.32e2_dp/0.15e2_dp*t142*t144 + REAL(2*t64*t147, KIND=dp) + t150 &
586 77657996 : *t118
587 77657996 : t155 = t29*t31*C
588 77657996 : t156 = t73*t116
589 : t159 = t126 + 2*D*E + t14*t140 + t115*t152 + t155*t156* &
590 77657996 : t118
591 77657996 : t162 = t35**2
592 77657996 : t163 = t162*t1
593 77657996 : t164 = t6*t5
594 77657996 : t167 = t61*t103*t104
595 77657996 : t170 = t34*E
596 77657996 : t173 = -0.16e2_dp/0.15e2_dp*t167*t109 + REAL(2*t64*t170, KIND=dp)
597 77657996 : t175 = t34*t103
598 77657996 : t176 = t64*t175
599 77657996 : t177 = t104*t106
600 77657996 : t178 = t177*t143
601 77657996 : t181 = E*t116
602 77657996 : t183 = -0.32e2_dp/0.15e2_dp*t176*t178 + t181*t118
603 77657996 : t187 = t104*t87*t119
604 : t190 = t61*E + t14*t173 + t115*t183 - 0.16e2_dp/0.15e2_dp*t115 &
605 77657996 : *t103*t187
606 : t194 = 2*E + t60 + t61*B + t14*t113 + t115*t123 + t125*t37 &
607 77657996 : *t159 + 3*t163*t164*t190
608 77657996 : t195 = t58*t194
609 77657996 : t196 = D*t35
610 77657996 : 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 77657996 : 0.9e1_dp*A*t54 - 0.4e1_dp/0.9e1_dp*t195*t199
613 77657996 : e_0 = e_0 + (t45*t202*Clda)*sx
614 : END IF
615 77657996 : IF (order >= 1 .OR. order >= -1) THEN
616 77657996 : t209 = rho**2
617 77657996 : srho = -t3/t164*t8*t4/0.18e2_dp - t3*t7/t209/0.6e1_dp
618 :
619 77657996 : t214 = t2*t7
620 77657996 : sndrho = t214*t8/0.6e1_dp
621 77657996 : t216 = t11*t13
622 77657996 : t217 = t216*t40
623 77657996 : t218 = dsdrho(rho, ndrho)
624 77657996 : t222 = 2*t217*t125*t37*t218
625 77657996 : t223 = a1*t11
626 77657996 : t224 = t13*t218
627 77657996 : t227 = t12*t11
628 77657996 : t228 = a2*t227
629 77657996 : t229 = t19*t218
630 77657996 : t232 = 2*t223*t224 + 4*t228*t229
631 77657996 : t234 = t14*t232*t39
632 77657996 : t235 = t21*t118
633 77657996 : t236 = t14*t235
634 77657996 : t237 = a3*t227
635 77657996 : t240 = a4*t17
636 77657996 : t244 = a5*t25
637 77657996 : t248 = 4*t237*t229 + 5*t240*t27*t218 + 6*t244*t31*t218
638 77657996 : t251 = t236*t125*t37*t248
639 77657996 : t255 = 0.2e1_dp/0.3e1_dp*t50*t125*t7*t4
640 77657996 : qrho = t222 + t234 - t251 + t255
641 77657996 : t256 = t216*t21
642 77657996 : t257 = t34*t41
643 77657996 : t261 = t232*t34
644 77657996 : t262 = t261*t41
645 77657996 : t265 = t118*t41
646 : prho = 0.9e1_dp/0.2e1_dp*t256*t257*t218 + 0.9e1_dp/0.4e1_dp*t14*t262 &
647 77657996 : - 0.9e1_dp/0.4e1_dp*t22*t265*t248
648 77657996 : t269 = dsdndrho(rho)
649 77657996 : t273 = t13*t269
650 77657996 : t276 = t19*t269
651 77657996 : t279 = 2*t223*t273 + 4*t228*t276
652 77657996 : t280 = t279*t34
653 77657996 : t281 = t280*t41
654 77657996 : 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 77657996 : t281 - 0.9e1_dp/0.4e1_dp*t22*t265*t292
657 : qndrho = 2*t217*t125*t37*t269 + t14*t279*t39 - t236*t125 &
658 77657996 : *t37*t292
659 77657996 : t308 = dexeirho(P, Q, Prho, Qrho)
660 77657996 : t312 = EXP(-q)
661 77657996 : t314 = t312*t106*t107
662 77657996 : t318 = 0.1e1_dp/t35
663 77657996 : t320 = 0.1e1_dp/t37
664 77657996 : t322 = 0.1e1_dp/t21*t33*t318*t1*t320
665 77657996 : t331 = 2*t216*t40*t218 + t14*t261 - t14*t235*t248
666 77657996 : t334 = t214*t4
667 77657996 : t339 = EXP(-t53)
668 77657996 : t341 = 0.1e1_dp/t51
669 77657996 : t347 = t56**2
670 77657996 : t349 = 0.1e1_dp/t347*t194
671 77657996 : t359 = D*t232
672 77657996 : t362 = t118*B
673 77657996 : t368 = t118*t248
674 77657996 : t370 = F1*t232*t34 - t71*t368
675 77657996 : t373 = t73*t11
676 77657996 : t382 = B*t51
677 77657996 : t385 = A*t56
678 77657996 : t393 = 0.1e1_dp/t86/t347
679 77657996 : t401 = t92*t93
680 77657996 : t402 = SQRT(0.31415926535897932385e1_dp)
681 77657996 : t404 = EXP(-p)
682 77657996 : t405 = 0.1e1_dp/t402*t404
683 77657996 : 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 77657996 : *(t262 - t235*t41*t248))
691 77657996 : t421 = t420*t104
692 77657996 : t424 = 0.1e1_dp/t227
693 77657996 : t425 = t105*t424
694 77657996 : t426 = t143*t218
695 77657996 : t429 = t86*t56
696 77657996 : t430 = t107*t429
697 77657996 : t431 = t430*t331
698 77657996 : t437 = t227*t19
699 77657996 : t445 = 0.1e1_dp/t117/t33
700 77657996 : t446 = t116*t445
701 77657996 : t451 = t121*t248
702 77657996 : t473 = t424*t107
703 77657996 : t475 = t473*t87*t218
704 77657996 : t479 = t108*t429*t331
705 77657996 : t484 = t118*C
706 77657996 : t497 = t105*t473
707 77657996 : t498 = t87*t21
708 77657996 : t503 = t105*t108
709 77657996 : t504 = t429*t21
710 77657996 : t517 = t64*t118
711 77657996 : t523 = C*t21
712 77657996 : t524 = t118*t232
713 77657996 : t527 = t445*t248
714 77657996 : t533 = t25*t31*C
715 77657996 : t534 = t118*t218
716 77657996 : t541 = t73*t21
717 77657996 : t568 = t118*E
718 77657996 : t581 = t64*t118*t103
719 77657996 : t590 = t104*t424
720 77657996 : t603 = t437*t105
721 77657996 : t604 = t87*t116
722 77657996 : t611 = t115*t105
723 77657996 : 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 77657996 : Clda)*sx
762 77657996 : t640 = dexeindrho(P, Q, Pndrho, Qndrho)
763 77657996 : t653 = 2*t216*t40*t269 + t14*t280 - t14*t235*t292
764 77657996 : t667 = D*t279
765 77657996 : t675 = t118*t292
766 77657996 : 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 77657996 : t409*(t281 - t235*t41*t292))
774 77657996 : t717 = t716*t104
775 77657996 : t720 = t143*t269
776 77657996 : t723 = t430*t653
777 77657996 : t739 = t121*t292
778 77657996 : t758 = t473*t87*t269
779 77657996 : t762 = t108*t429*t653
780 77657996 : t800 = t118*t279
781 77657996 : t803 = t445*t292
782 77657996 : 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 77657996 : + 0.4e1_dp/0.9e1_dp*t195*qndrho*t199)*Clda)*sx
817 : END IF
818 :
819 77657996 : 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 32423 : 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 32423 : IF (order >= 0) THEN
856 32423 : t1 = 3**(0.1e1_dp/0.3e1_dp)
857 32423 : t2 = t1**2
858 32423 : t3 = ndrho*t2
859 32423 : t4 = pi**2
860 32423 : t5 = t4*rho
861 32423 : t6 = t5**(0.1e1_dp/0.3e1_dp)
862 32423 : t7 = 0.1e1_dp/t6
863 32423 : t8 = 0.1e1_dp/rho
864 32423 : ssval = t3*t7*t8/0.6e1_dp
865 :
866 32423 : t11 = red(rho, ndrho)
867 32423 : t12 = t11**2
868 32423 : t13 = sscale**2
869 32423 : t14 = t12*t13
870 32423 : t17 = t12**2
871 32423 : t19 = t13**2
872 32423 : t21 = a1*t12*t13 + a2*t17*t19
873 32423 : t22 = t14*t21
874 32423 : t25 = t17*t11
875 32423 : t27 = t19*sscale
876 32423 : t29 = t17*t12
877 32423 : t31 = t19*t13
878 32423 : t33 = 1 + a3*t17*t19 + a4*t25*t27 + a5*t29*t31
879 32423 : t34 = 0.1e1_dp/t33
880 32423 : t35 = R**2
881 32423 : t37 = t6**2
882 32423 : t38 = t2*t37
883 32423 : t39 = t34*t35*t38
884 32423 : q = t22*t39
885 32423 : t40 = t21*t34
886 32423 : t41 = 0.1e1_dp/A
887 32423 : p = 0.9e1_dp/0.4e1_dp*t14*t40*t41
888 32423 : t44 = rho**(0.1e1_dp/0.3e1_dp)
889 32423 : t45 = t44*rho
890 32423 : t46 = exei(P, Q)
891 32423 : t48 = expint(1, q)
892 32423 : t50 = t14*t40
893 32423 : t51 = D + t50
894 32423 : t52 = t51*t35
895 32423 : t53 = t52*t38
896 32423 : t54 = expint(1, t53)
897 32423 : t56 = t51**2
898 32423 : t58 = 0.1e1_dp/t56/t51
899 32423 : t60 = D*C
900 32423 : t61 = D**2
901 32423 : t64 = D*t21
902 32423 : t65 = t34*B
903 32423 : t74 = g1 + g2*t12*t13 + g3*t17*t19 + g4*t25*t27
904 32423 : t75 = E*t74
905 32423 : t77 = F1*t21
906 32423 : t79 = t77*t34 + F2
907 32423 : t81 = t40*C + 2*t64*t65 + 2*t75 + t60*t79
908 32423 : t83 = t17*t19
909 32423 : t84 = t21**2
910 32423 : t85 = t33**2
911 32423 : t86 = 0.1e1_dp/t85
912 32423 : t89 = C*t79
913 32423 : t91 = t84*t86*B + t40*t89
914 32423 : t93 = t35*t2
915 32423 : t94 = t61*C
916 32423 : t95 = D*E
917 32423 : t97 = E*t21
918 32423 : t102 = t34*C
919 32423 : t106 = 2*t97*t34 + 2*t95*t74 + 2*t64*t102 + t94*t79
920 32423 : t110 = t102*t79
921 32423 : t113 = C*t84
922 32423 : t115 = 2*t75*t40 + 2*t64*t110 + t113*t86
923 32423 : t117 = t29*t31
924 32423 : t118 = t117*C
925 32423 : t119 = t79*t84
926 32423 : t122 = t94 + 2*t95 + t14*t106 + t83*t115 + t118*t119*t86
927 32423 : t125 = t35**2
928 32423 : t126 = t125*t1
929 32423 : t127 = t6*t5
930 32423 : t128 = t61*E
931 32423 : t130 = t34*E
932 32423 : t133 = t128*t74 + 2*t64*t130
933 32423 : t135 = t130*t74
934 32423 : t138 = E*t84
935 32423 : t140 = 2*t64*t135 + t138*t86
936 32423 : t142 = t117*E
937 32423 : t143 = t74*t84
938 32423 : t146 = t128 + t14*t133 + t83*t140 + t142*t143*t86
939 : t150 = 2*E + t60 + t61*B + t14*t81 + t83*t91 + t93*t37* &
940 32423 : t122 + 3*t126*t127*t146
941 32423 : t151 = t58*t150
942 32423 : t152 = D*t35
943 32423 : 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 32423 : 0.9e1_dp*A*t54 - 0.4e1_dp/0.9e1_dp*t151*t155
946 32423 : e_0 = e_0 + (t45*t158*Clda)*sx
947 : END IF
948 32423 : IF (order >= 1 .OR. order == -1) THEN
949 27653 : t165 = rho**2
950 27653 : srho = -t3/t127*t8*t4/0.18e2_dp - t3*t7/t165/0.6e1_dp
951 27653 : t170 = t2*t7
952 27653 : sndrho = t170*t8/0.6e1_dp
953 27653 : t172 = t11*t13
954 27653 : t173 = t172*t40
955 27653 : t174 = dsdrho(rho, ndrho)
956 27653 : t178 = 2*t173*t93*t37*t174
957 27653 : t179 = a1*t11
958 27653 : t180 = t13*t174
959 27653 : t183 = t12*t11
960 27653 : t184 = a2*t183
961 27653 : t185 = t19*t174
962 27653 : t188 = 2*t179*t180 + 4*t184*t185
963 27653 : t190 = t14*t188*t39
964 27653 : t191 = t21*t86
965 27653 : t192 = t14*t191
966 27653 : t193 = a3*t183
967 27653 : t196 = a4*t17
968 27653 : t197 = t27*t174
969 27653 : t200 = a5*t25
970 27653 : t204 = 4*t193*t185 + 5*t196*t197 + 6*t200*t31*t174
971 27653 : t207 = t192*t93*t37*t204
972 27653 : t211 = 0.2e1_dp/0.3e1_dp*t50*t93*t7*t4
973 27653 : qrho = t178 + t190 - t207 + t211
974 27653 : t212 = t172*t21
975 27653 : t213 = t34*t41
976 27653 : t217 = t188*t34
977 27653 : t221 = t86*t41
978 : prho = 0.9e1_dp/0.2e1_dp*t212*t213*t174 + 0.9e1_dp/0.4e1_dp*t14*t217 &
979 27653 : *t41 - 0.9e1_dp/0.4e1_dp*t22*t221*t204
980 27653 : t225 = dsdndrho(rho)
981 27653 : t229 = t13*t225
982 27653 : t232 = t19*t225
983 27653 : t235 = 2*t179*t229 + 4*t184*t232
984 27653 : t236 = t235*t34
985 27653 : t242 = t27*t225
986 27653 : 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 27653 : t236*t41 - 0.9e1_dp/0.4e1_dp*t22*t221*t248
989 : qndrho = 2*t173*t93*t37*t225 + t14*t235*t39 - t192*t93 &
990 27653 : *t37*t248
991 27653 : t264 = dexeirho(P, Q, Prho, Qrho)
992 27653 : t268 = EXP(-q)
993 27653 : t272 = t268/t12/t13
994 27653 : t276 = 0.1e1_dp/t35
995 27653 : t278 = 0.1e1_dp/t37
996 27653 : t280 = 0.1e1_dp/t21*t33*t276*t1*t278
997 27653 : t287 = t191*t204
998 27653 : t289 = 2*t172*t40*t174 + t14*t217 - t14*t287
999 27653 : t292 = t170*t4
1000 27653 : t297 = EXP(-t53)
1001 27653 : t299 = 0.1e1_dp/t51
1002 27653 : t305 = t56**2
1003 27653 : t307 = 0.1e1_dp/t305*t150
1004 27653 : t317 = D*t188
1005 27653 : t320 = t86*B
1006 27653 : t324 = g2*t11
1007 27653 : t327 = g3*t183
1008 27653 : t330 = g4*t17
1009 27653 : t333 = 2*t324*t180 + 4*t327*t185 + 5*t330*t197
1010 27653 : t334 = E*t333
1011 27653 : t338 = t86*t204
1012 27653 : t340 = F1*t188*t34 - t77*t338
1013 27653 : t344 = t183*t19
1014 27653 : t352 = 0.1e1_dp/t85/t33
1015 27653 : t353 = t84*t352
1016 27653 : t358 = t89*t204
1017 27653 : t380 = t86*C
1018 27653 : t394 = t64*t86
1019 27653 : t398 = C*t21
1020 27653 : t399 = t86*t188
1021 27653 : t401 = t352*t204
1022 27653 : t406 = t25*t31
1023 27653 : t407 = t406*C
1024 27653 : t408 = t86*t174
1025 27653 : t415 = t79*t21
1026 27653 : t435 = t86*E
1027 27653 : t454 = t406*E
1028 27653 : 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 27653 : Clda)*sx
1053 27653 : t485 = dexeindrho(P, Q, Pndrho, Qndrho)
1054 27653 : t496 = t191*t248
1055 27653 : t498 = 2*t172*t40*t225 + t14*t236 - t14*t496
1056 27653 : t512 = D*t235
1057 27653 : t524 = 2*t324*t229 + 4*t327*t232 + 5*t330*t242
1058 27653 : t525 = E*t524
1059 27653 : t529 = t86*t248
1060 27653 : t531 = F1*t235*t34 - t77*t529
1061 27653 : t545 = t89*t248
1062 27653 : t579 = t86*t235
1063 27653 : t581 = t352*t248
1064 27653 : 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 27653 : *t155)*Clda)*sx
1085 : END IF
1086 :
1087 32423 : 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 233061717 : 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 233061717 : exei = 0.0_dp
1112 233061717 : IF (P < expcutoff) THEN
1113 : !Use exact product
1114 233061717 : IF (P + Q < 0.5_dp) THEN
1115 11768955 : tmp = -euler - LOG(P + Q) + P + Q
1116 11768955 : 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 11768955 : tmp = tmp + 1.0_dp/600.0_dp*(P + Q)**5
1118 11768955 : exei = EXP(P)*tmp
1119 : ELSE
1120 221292762 : exei = EXP(P)*expint(1, Q + P)
1121 : END IF
1122 : ELSE
1123 : !Use approximation
1124 0 : tmp = 1.0_dp/P
1125 : ! *** 1st order
1126 0 : exei = tmp
1127 : ! *** 2nd order
1128 0 : tmp = tmp/P
1129 0 : exei = exei - (Q + 1.0_dp)*tmp
1130 : ! *** 3rd order
1131 0 : tmp = tmp/P
1132 0 : Q2 = Q*Q
1133 0 : exei = exei + (2.0_dp*Q + Q2 + 2.0_dp)*tmp
1134 : ! *** 4th order
1135 0 : tmp = tmp/P
1136 0 : Q3 = Q2*Q
1137 0 : exei = exei - (3.0_dp*Q2 + 6.0_dp*Q + Q3 + 6.0_dp)*tmp
1138 : ! *** 5th order
1139 0 : tmp = tmp/P
1140 0 : Q4 = Q3*Q
1141 0 : exei = exei + (24.0_dp + 4.0_dp*Q3 + Q4 + 12.0_dp*Q2 + 24.0_dp*Q)*tmp
1142 :
1143 : ! *** scaling factor
1144 0 : exei = EXP(-Q)*exei
1145 : END IF
1146 233061717 : END FUNCTION exei
1147 :
1148 : ! **************************************************************************************************
1149 : !> \brief ...
1150 : !> \param P ...
1151 : !> \param Q ...
1152 : !> \param dPrho ...
1153 : !> \param dQrho ...
1154 : !> \return ...
1155 : ! **************************************************************************************************
1156 77685649 : ELEMENTAL FUNCTION dexeirho(P, Q, dPrho, dQrho)
1157 : REAL(dp), INTENT(IN) :: P, Q, dPrho, dQrho
1158 : REAL(dp) :: dexeirho
1159 :
1160 77685649 : dexeirho = dPrho*exei(P, Q) - (dPrho + dQrho)/(P + Q)*EXP(-Q)
1161 77685649 : END FUNCTION dexeirho
1162 :
1163 : ! **************************************************************************************************
1164 : !> \brief ...
1165 : !> \param P ...
1166 : !> \param Q ...
1167 : !> \param dPndrho ...
1168 : !> \param dQndrho ...
1169 : !> \return ...
1170 : ! **************************************************************************************************
1171 77685649 : ELEMENTAL FUNCTION dexeindrho(P, Q, dPndrho, dQndrho)
1172 : REAL(dp), INTENT(IN) :: P, Q, dPndrho, dQndrho
1173 : REAL(dp) :: dexeindrho
1174 :
1175 77685649 : dexeindrho = dPndrho*exei(P, Q) - (dPndrho + dQndrho)/(P + Q)*EXP(-Q)
1176 77685649 : END FUNCTION dexeindrho
1177 :
1178 : ! **************************************************************************************************
1179 : !> \brief ...
1180 : !> \param rho ...
1181 : !> \param ndrho ...
1182 : !> \return ...
1183 : ! **************************************************************************************************
1184 77690419 : ELEMENTAL FUNCTION red(rho, ndrho)
1185 : REAL(dp), INTENT(IN) :: rho, ndrho
1186 : REAL(dp) :: red
1187 :
1188 77690419 : red = 1.0_dp/6.0_dp*ndrho*3.0_dp**(2.0_dp/3.0_dp)
1189 77690419 : red = red/(pi**(2.0_dp/3.0_dp))
1190 77690419 : red = red*MAX(1.0_dp/rho**(4.0_dp/3.0_dp), EPS_RHO)
1191 77690419 : END FUNCTION red
1192 :
1193 : ! **************************************************************************************************
1194 : !> \brief ...
1195 : !> \param rho ...
1196 : !> \param ndrho ...
1197 : !> \return ...
1198 : ! **************************************************************************************************
1199 77685649 : ELEMENTAL FUNCTION dsdrho(rho, ndrho)
1200 : REAL(dp), INTENT(IN) :: rho, ndrho
1201 : REAL(dp) :: dsdrho
1202 :
1203 77685649 : dsdrho = -2.0_dp/9.0_dp*ndrho*3.0**(2.0_dp/3.0_dp)
1204 77685649 : dsdrho = dsdrho/(pi**(2.0_dp/3.0_dp))
1205 77685649 : dsdrho = dsdrho*MAX(1.0_dp/rho**(7.0_dp/3.0_dp), EPS_RHO)
1206 77685649 : END FUNCTION dsdrho
1207 :
1208 : ! **************************************************************************************************
1209 : !> \brief ...
1210 : !> \param rho ...
1211 : !> \return ...
1212 : ! **************************************************************************************************
1213 77685649 : 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 77685649 : dsdndrho = dsdndrho/(pi**(2.0_dp/3.0_dp))
1219 77685649 : dsdndrho = dsdndrho*MAX(1.0_dp/rho**(4.0_dp/3.0_dp), EPS_RHO)
1220 77685649 : 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 376673600 : 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 376673600 : nm1 = n - 1
1246 :
1247 376673600 : 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 376673600 : ELSE IF (n .EQ. 0) THEN !Special case.
1250 0 : expint = EXP(-x)/x
1251 376673600 : ELSE IF (x .EQ. 0.0_dp) THEN !Another special case.
1252 0 : expint = 1.0_dp/nm1
1253 376673600 : ELSE IF (x .GT. 1.0_dp) THEN !Lentz's algorithm (5.2).
1254 244309972 : b = x + n
1255 244309972 : c = 1.0_dp/FPMIN
1256 244309972 : d = 1.0_dp/b
1257 244309972 : h = d
1258 7648239479 : DO i = 1, MAXIT
1259 7648239479 : a = -i*(nm1 + i)
1260 7648239479 : b = b + 2.0_dp
1261 7648239479 : d = 1.0_dp/(a*d + b)
1262 7648239479 : c = b + a/c
1263 7648239479 : del = c*d
1264 7648239479 : h = h*del
1265 7648239479 : IF (ABS(del - 1.0_dp) .LT. EPS .OR. i == MAXIT) THEN
1266 244309972 : expint = h*EXP(-x)
1267 244309972 : RETURN
1268 : END IF
1269 : END DO
1270 : ELSE !Evaluate series.
1271 132363628 : IF (nm1 .NE. 0) THEN !Set first term.
1272 0 : expint = 1.0_dp/nm1
1273 : ELSE
1274 132363628 : expint = -LOG(x) - euler
1275 : END IF
1276 : fact = 1.0_dp
1277 1136765152 : DO i = 1, MAXIT
1278 1136765152 : fact = -fact*x/i
1279 1136765152 : IF (i .NE. nm1) THEN
1280 1136765152 : del = -fact/(i - nm1)
1281 : ELSE
1282 : psi = -euler !Compute I(n).
1283 0 : DO ii = 1, nm1
1284 0 : psi = psi + 1.0_dp/ii
1285 : END DO
1286 0 : del = fact*(-LOG(x) + psi)
1287 : END IF
1288 1136765152 : expint = expint + del
1289 1136765152 : 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 :
|
