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