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 Calculate the LDA functional according to Vosk, Wilk and Nusair
10 : !> Literature: S. H. Vosko, L. Wilk and M. Nusair,
11 : !> Can. J. Phys. 58, 1200 (1980)
12 : !> \note
13 : !> Order of derivatives is: LDA 0; 1; 2; 3;
14 : !> \par History
15 : !> JGH (26.02.2003) : OpenMP enabled
16 : !> fawzi (04.2004) : adapted to the new xc interface
17 : !> MG 01.2007 : added scaling
18 : !> MG 01.2007 : bug fix in vwn_lda_xx
19 : !> \author JGH (26.02.2002)
20 : ! **************************************************************************************************
21 : MODULE xc_vwn
22 : USE bibliography, ONLY: Vosko1980,&
23 : cite_reference
24 : USE input_section_types, ONLY: section_vals_type,&
25 : section_vals_val_get
26 : USE kinds, ONLY: dp
27 : USE xc_derivative_desc, ONLY: 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_functionals_utilities, ONLY: calc_srs_pw,&
35 : set_util
36 : USE xc_input_constants, ONLY: do_vwn3,&
37 : do_vwn5
38 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
39 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
40 : xc_rho_set_type
41 : #include "../base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 :
45 : PRIVATE
46 :
47 : ! *** Global parameters ***
48 :
49 : REAL(KIND=dp), PARAMETER :: a = 0.0310907_dp, &
50 : af = 0.01554535_dp, &
51 : aa = -0.0168868639404_dp !-1/(6pi^2)
52 :
53 : PUBLIC :: vwn_lda_info, vwn_lda_eval, vwn_lsd_eval, vwn_lsd_info
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_vwn'
56 :
57 : REAL(KIND=dp) :: eps_rho, b, c, x0, bf, cf, x0f, ba, ca, x0a
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief ...
63 : !> \param cutoff ...
64 : !> \param vwn_params ...
65 : ! **************************************************************************************************
66 1163 : SUBROUTINE vwn_init(cutoff, vwn_params)
67 :
68 : REAL(KIND=dp), INTENT(IN) :: cutoff
69 : TYPE(section_vals_type), POINTER :: vwn_params
70 :
71 : INTEGER :: TYPE
72 :
73 1163 : CALL section_vals_val_get(vwn_params, "functional_type", i_val=TYPE)
74 1163 : eps_rho = cutoff
75 1163 : CALL set_util(cutoff)
76 :
77 1163 : CALL cite_reference(Vosko1980)
78 :
79 1109 : SELECT CASE (TYPE)
80 : CASE (do_vwn5)
81 1109 : b = 3.72744_dp
82 1109 : c = 12.9352_dp
83 1109 : x0 = -0.10498_dp
84 1109 : bf = 7.06042_dp
85 1109 : cf = 18.0578_dp
86 1109 : x0f = -0.32500_dp
87 1109 : ba = 1.13107_dp
88 1109 : ca = 13.0045_dp
89 1109 : x0a = -0.00475840_dp
90 : CASE (do_vwn3)
91 54 : b = 13.0720_dp
92 54 : c = 42.7198_dp
93 54 : x0 = -0.409286_dp
94 54 : bf = 20.1231_dp
95 54 : cf = 101.578_dp
96 54 : x0f = -0.743294_dp
97 54 : ba = 1.13107_dp
98 54 : ca = 13.0045_dp
99 54 : x0a = -0.00475840_dp
100 : CASE DEFAULT
101 1163 : CPABORT(" Only functionals VWN3 and VWN5 are supported")
102 :
103 : END SELECT
104 :
105 1163 : END SUBROUTINE vwn_init
106 :
107 : ! **************************************************************************************************
108 : !> \brief ...
109 : !> \param reference ...
110 : !> \param shortform ...
111 : !> \param needs ...
112 : !> \param max_deriv ...
113 : ! **************************************************************************************************
114 946 : SUBROUTINE vwn_lda_info(reference, shortform, needs, max_deriv)
115 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
116 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
117 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
118 :
119 946 : IF (PRESENT(reference)) THEN
120 : reference = "S. H. Vosko, L. Wilk and M. Nusair,"// &
121 23 : " Can. J. Phys. 58, 1200 (1980) {LDA version}"
122 : END IF
123 946 : IF (PRESENT(shortform)) THEN
124 23 : shortform = "Vosko-Wilk-Nusair Functional {LDA}"
125 : END IF
126 946 : IF (PRESENT(needs)) THEN
127 923 : needs%rho = .TRUE.
128 : END IF
129 946 : IF (PRESENT(max_deriv)) max_deriv = 3
130 :
131 946 : END SUBROUTINE vwn_lda_info
132 :
133 : ! **************************************************************************************************
134 : !> \brief ...
135 : !> \param reference ...
136 : !> \param shortform ...
137 : !> \param needs ...
138 : !> \param max_deriv ...
139 : ! **************************************************************************************************
140 20 : SUBROUTINE vwn_lsd_info(reference, shortform, needs, max_deriv)
141 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
142 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
143 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
144 :
145 20 : IF (PRESENT(reference)) THEN
146 : reference = "S. H. Vosko, L. Wilk and M. Nusair,"// &
147 2 : " Can. J. Phys. 58, 1200 (1980) {LDA version}"
148 : END IF
149 20 : IF (PRESENT(shortform)) THEN
150 2 : shortform = "Vosko-Wilk-Nusair Functional {LDA}"
151 : END IF
152 20 : IF (PRESENT(needs)) THEN
153 18 : needs%rho_spin = .TRUE.
154 : END IF
155 20 : IF (PRESENT(max_deriv)) max_deriv = 3
156 :
157 20 : END SUBROUTINE vwn_lsd_info
158 :
159 : ! **************************************************************************************************
160 : !> \brief ...
161 : !> \param rho_set ...
162 : !> \param deriv_set ...
163 : !> \param order ...
164 : !> \param vwn_params ...
165 : ! **************************************************************************************************
166 1131 : SUBROUTINE vwn_lda_eval(rho_set, deriv_set, order, vwn_params)
167 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
168 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
169 : INTEGER, INTENT(in) :: order
170 : TYPE(section_vals_type), POINTER :: vwn_params
171 :
172 : CHARACTER(len=*), PARAMETER :: routineN = 'vwn_lda_eval'
173 :
174 : INTEGER :: handle, npoints
175 : INTEGER, DIMENSION(2, 3) :: bo
176 : REAL(KIND=dp) :: epsilon_rho, sc
177 1131 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x
178 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
179 1131 : POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, rho
180 : TYPE(xc_derivative_type), POINTER :: deriv
181 :
182 1131 : CALL timeset(routineN, handle)
183 :
184 1131 : CALL section_vals_val_get(vwn_params, "scale_c", r_val=sc)
185 :
186 : CALL xc_rho_set_get(rho_set, rho=rho, &
187 1131 : local_bounds=bo, rho_cutoff=epsilon_rho)
188 1131 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
189 1131 : CALL vwn_init(epsilon_rho, vwn_params)
190 :
191 3393 : ALLOCATE (x(npoints))
192 1131 : CALL calc_srs_pw(rho, x, npoints)
193 :
194 1131 : IF (order >= 1) THEN
195 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
196 1129 : allocate_deriv=.TRUE.)
197 1129 : CALL xc_derivative_get(deriv, deriv_data=e_0)
198 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
199 1129 : allocate_deriv=.TRUE.)
200 1129 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
201 :
202 1129 : CALL vwn_lda_01(rho, x, e_0, e_rho, npoints, sc)
203 : ELSE
204 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
205 2 : allocate_deriv=.TRUE.)
206 2 : CALL xc_derivative_get(deriv, deriv_data=e_0)
207 :
208 2 : CALL vwn_lda_0(rho, x, e_0, npoints, sc)
209 : END IF
210 1131 : IF (order >= 2) THEN
211 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
212 0 : allocate_deriv=.TRUE.)
213 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
214 :
215 0 : CALL vwn_lda_2(rho, x, e_rho_rho, npoints, sc)
216 : END IF
217 1131 : IF (order >= 3) THEN
218 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
219 0 : allocate_deriv=.TRUE.)
220 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
221 :
222 0 : CALL vwn_lda_3(rho, x, e_rho_rho_rho, npoints, sc)
223 : END IF
224 1131 : IF (order > 3) THEN
225 0 : CPABORT("derivatives bigger than 3 not implemented")
226 : END IF
227 :
228 1131 : DEALLOCATE (x)
229 1131 : CALL timestop(handle)
230 :
231 1131 : END SUBROUTINE vwn_lda_eval
232 :
233 : ! **************************************************************************************************
234 : !> \brief ...
235 : !> \param rho ...
236 : !> \param x ...
237 : !> \param e_0 ...
238 : !> \param npoints ...
239 : !> \param sc ...
240 : ! **************************************************************************************************
241 2 : SUBROUTINE vwn_lda_0(rho, x, e_0, npoints, sc)
242 :
243 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, x
244 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
245 : INTEGER, INTENT(in) :: npoints
246 : REAL(KIND=dp) :: sc
247 :
248 : INTEGER :: ip
249 : REAL(KIND=dp) :: at, dpx, ln1, ln2, px, px0, q, xb
250 :
251 2 : q = SQRT(4.0_dp*c - b*b)
252 2 : xb = 2.0_dp*x0 + b
253 2 : px0 = x0*x0 + b*x0 + c
254 :
255 : !$OMP PARALLEL DO DEFAULT(NONE) &
256 : !$OMP SHARED (npoints, rho, eps_rho, x, c, b, e_0, sc) &
257 : !$OMP SHARED (q, x0, px0, xb) &
258 2 : !$OMP PRIVATE(ip,px,dpx,at,ln1,ln2)
259 :
260 : DO ip = 1, npoints
261 :
262 : IF (rho(ip) > eps_rho) THEN
263 : px = x(ip)*x(ip) + b*x(ip) + c
264 : dpx = 2.0_dp*x(ip) + b
265 : at = 2.0_dp/q*ATAN(q/dpx)
266 : ln1 = LOG(x(ip)*x(ip)/px)
267 : ln2 = LOG((x(ip) - x0)**2/px)
268 : e_0(ip) = e_0(ip) + a*(ln1 + b*at - b*x0/px0*(ln2 + xb*at))*rho(ip)*sc
269 : END IF
270 :
271 : END DO
272 :
273 : !$OMP END PARALLEL DO
274 :
275 2 : END SUBROUTINE vwn_lda_0
276 :
277 : ! **************************************************************************************************
278 : !> \brief ...
279 : !> \param rho ...
280 : !> \param x ...
281 : !> \param e_0 ...
282 : !> \param e_rho ...
283 : !> \param npoints ...
284 : !> \param sc ...
285 : ! **************************************************************************************************
286 1129 : SUBROUTINE vwn_lda_01(rho, x, e_0, e_rho, npoints, sc)
287 :
288 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, x
289 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho
290 : INTEGER, INTENT(in) :: npoints
291 : REAL(KIND=dp) :: sc
292 :
293 : INTEGER :: ip
294 : REAL(KIND=dp) :: at, dat, dex, dln1, dln2, dpx, ex, ln1, &
295 : ln2, pa, px, px0, q, xb
296 :
297 1129 : q = SQRT(4.0_dp*c - b*b)
298 1129 : xb = 2.0_dp*x0 + b
299 1129 : px0 = x0*x0 + b*x0 + c
300 :
301 : !$OMP PARALLEL DO DEFAULT(NONE) &
302 : !$OMP SHARED(npoints, rho, eps_rho, x, b, c, e_0, sc) &
303 : !$OMP SHARED(x0, q, xb, px0, e_rho) &
304 1129 : !$OMP PRIVATE(ip,px,dpx,at,pa,dat,ln1,dln1,ln2,dln2,ex,dex)
305 :
306 : DO ip = 1, npoints
307 :
308 : IF (rho(ip) > eps_rho) THEN
309 : px = x(ip)*x(ip) + b*x(ip) + c
310 : dpx = 2.0_dp*x(ip) + b
311 : at = 2.0_dp/q*ATAN(q/dpx)
312 : pa = 4.0_dp*x(ip)*x(ip) + 4.0_dp*b*x(ip) + b*b + q*q
313 : dat = -4.0_dp/pa
314 : ln1 = LOG(x(ip)*x(ip)/px)
315 : dln1 = (b*x(ip) + 2.0_dp*c)/(x(ip)*px)
316 : ln2 = LOG((x(ip) - x0)**2/px)
317 : dln2 = (x(ip)*(b + 2.0_dp*x0) + 2.0_dp*c + x0*b)/((x(ip) - x0)*px)
318 : ex = a*(ln1 + b*at - b*x0/px0*(ln2 + xb*at))
319 : dex = a*(dln1 + b*dat - b*x0/px0*(dln2 + xb*dat))
320 : e_0(ip) = e_0(ip) + ex*rho(ip)*sc
321 : e_rho(ip) = e_rho(ip) + (ex - x(ip)*dex/6.0_dp)*sc
322 : END IF
323 :
324 : END DO
325 :
326 : !$OMP END PARALLEL DO
327 :
328 1129 : END SUBROUTINE vwn_lda_01
329 :
330 : ! **************************************************************************************************
331 : !> \brief ...
332 : !> \param rho ...
333 : !> \param x ...
334 : !> \param e_rho_rho ...
335 : !> \param npoints ...
336 : !> \param sc ...
337 : ! **************************************************************************************************
338 0 : SUBROUTINE vwn_lda_2(rho, x, e_rho_rho, npoints, sc)
339 :
340 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, x
341 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
342 : INTEGER, INTENT(in) :: npoints
343 : REAL(KIND=dp) :: sc
344 :
345 : INTEGER :: ip
346 : REAL(KIND=dp) :: at, d2at, d2ex, d2ln1, d2ln2, dat, dex, &
347 : dln1, dln2, dpx, fp, ln1, ln2, pa, px, &
348 : px0, q, xb
349 :
350 0 : q = SQRT(4.0_dp*c - b*b)
351 0 : xb = 2.0_dp*x0 + b
352 0 : px0 = x0*x0 + b*x0 + c
353 0 : fp = -b*x0/px0
354 :
355 : !$OMP PARALLEL DO DEFAULT(NONE) &
356 : !$OMP PRIVATE(ip,px,dpx,at,pa,dat,d2at,ln1,dln1,d2ln1) &
357 : !$OMP PRIVATE(ln2,dln2,d2ln2,dex,d2ex) &
358 : !$OMP SHARED(npoints, rho, fp, eps_rho, x, b, c, q, x0) &
359 0 : !$OMP SHARED(xb, e_rho_rho, sc)
360 :
361 : DO ip = 1, npoints
362 :
363 : IF (rho(ip) > eps_rho) THEN
364 : px = x(ip)*x(ip) + b*x(ip) + c
365 : dpx = 2.0_dp*x(ip) + b
366 : at = 2.0_dp/q*ATAN(q/dpx)
367 : pa = 4.0_dp*x(ip)*x(ip) + 4.0_dp*b*x(ip) + b*b + q*q
368 : dat = -4.0_dp/pa
369 : d2at = 16.0_dp*dpx/(pa*pa)
370 : ln1 = LOG(x(ip)*x(ip)/px)
371 : dln1 = (b*x(ip) + 2.0_dp*c)/(x(ip)*px)
372 : d2ln1 = b/(x(ip)*px) - (b*x(ip) + 2.0_dp*c)/(x(ip)*px)**2*(px + x(ip)*dpx)
373 : ln2 = LOG((x(ip) - x0)**2/px)
374 : dln2 = (x(ip)*xb + 2.0_dp*c + x0*b)/((x(ip) - x0)*px)
375 : d2ln2 = xb/((x(ip) - x0)*px) - (x(ip)*xb + 2.0_dp*c + x0*b)/((x(ip) - x0)*px)**2 &
376 : *(px + (x(ip) - x0)*dpx)
377 : dex = a*(dln1 + b*dat + fp*(dln2 + xb*dat))
378 : d2ex = a*(d2ln1 + b*d2at + fp*(d2ln2 + xb*d2at))
379 : e_rho_rho(ip) = e_rho_rho(ip) &
380 : + (x(ip)/(36.0_dp*rho(ip))*(x(ip)*d2ex - 5.0_dp*dex))*sc
381 : END IF
382 :
383 : END DO
384 : !$OMP END PARALLEL DO
385 :
386 0 : END SUBROUTINE vwn_lda_2
387 :
388 : ! **************************************************************************************************
389 : !> \brief ...
390 : !> \param rho ...
391 : !> \param x ...
392 : !> \param e_rho_rho_rho ...
393 : !> \param npoints ...
394 : !> \param sc ...
395 : ! **************************************************************************************************
396 0 : SUBROUTINE vwn_lda_3(rho, x, e_rho_rho_rho, npoints, sc)
397 :
398 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, x
399 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
400 : INTEGER, INTENT(in) :: npoints
401 : REAL(KIND=dp) :: sc
402 :
403 : INTEGER :: ip
404 : REAL(KIND=dp) :: at, ax, bx, cx, d2at, d2bx, d2dx, d2ex, &
405 : d2ln1, d2ln2, d3at, d3ex, d3ln1, &
406 : d3ln2, dat, dbx, ddx, dex, dln1, dln2, &
407 : dpx, dx, fp, ln1, ln2, pa, px, px0, q, &
408 : xb
409 :
410 0 : q = SQRT(4.0_dp*c - b*b)
411 0 : xb = 2.0_dp*x0 + b
412 0 : px0 = x0*x0 + b*x0 + c
413 0 : fp = -b*x0/px0
414 :
415 : !$OMP PARALLEL DO DEFAULT(NONE) &
416 : !$OMP SHARED(npoints, rho, eps_rho, b, c, x, x0, sc) &
417 : !$OMP SHARED(q, xb, px0, fp, e_rho_rho_rho) &
418 : !$OMP PRIVATE(ip,px,dpx,at,pa,dat,d2at,d3at,ln1,ax,bx) &
419 : !$OMP PRIVATE(dbx,d2bx,dln1,d2ln1,d3ln1,ln2,cx,dx,ddx,d2dx) &
420 0 : !$OMP PRIVATE(dln2,d2ln2,d3ln2,dex,d2ex,d3ex)
421 : DO ip = 1, npoints
422 :
423 : IF (rho(ip) > eps_rho) THEN
424 : px = x(ip)*x(ip) + b*x(ip) + c
425 : dpx = 2.0_dp*x(ip) + b
426 : at = 2.0_dp/q*ATAN(q/dpx)
427 : pa = 4.0_dp*x(ip)*x(ip) + 4.0_dp*b*x(ip) + b*b + q*q
428 : dat = -4.0_dp/pa
429 : d2at = 16.0_dp*dpx/(pa*pa)
430 : d3at = 32.0_dp/(pa*pa)*(1.0_dp - 4.0_dp*dpx*dpx/pa)
431 : ln1 = LOG(x(ip)*x(ip)/px)
432 : ax = b*x(ip) + 2.0_dp*c
433 : bx = x(ip)*px
434 : dbx = px + x(ip)*dpx
435 : d2bx = 2.0_dp*(dpx + x(ip))
436 : dln1 = ax/bx
437 : d2ln1 = (b*bx - ax*dbx)/(bx*bx)
438 : d3ln1 = -ax*d2bx/(bx*bx) - 2.0_dp*d2ln1*dbx/bx
439 : ln2 = LOG((x(ip) - x0)**2/px)
440 : cx = x(ip)*xb + 2.0_dp*c + x0*b
441 : dx = (x(ip) - x0)*px
442 : ddx = px + (x(ip) - x0)*dpx
443 : d2dx = 2.0_dp*(dpx + (x(ip) - x0))
444 : dln2 = cx/dx
445 : d2ln2 = (xb*dx - cx*ddx)/(dx*dx)
446 : d3ln2 = -cx*d2dx/(dx*dx) - 2.0_dp*d2ln2*ddx/dx
447 : dex = a*(dln1 + b*dat + fp*(dln2 + xb*dat))
448 : d2ex = a*(d2ln1 + b*d2at + fp*(d2ln2 + xb*d2at))
449 : d3ex = a*(d3ln1 + b*d3at + fp*(d3ln2 + xb*d3at))
450 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
451 : - (7.0_dp*x(ip)/(216.0_dp*rho(ip)*rho(ip))*(x(ip)*d2ex - 5.0_dp*dex) + &
452 : x(ip)*x(ip)/(216.0_dp*rho(ip)*rho(ip))*(x(ip)*d3ex - 4.0_dp*d2ex))*sc
453 : END IF
454 :
455 : END DO
456 : !$OMP END PARALLEL DO
457 :
458 0 : END SUBROUTINE vwn_lda_3
459 :
460 : ! **************************************************************************************************
461 : !> \brief ...
462 : !> \param rho_set ...
463 : !> \param deriv_set ...
464 : !> \param order ...
465 : !> \param vwn_params ...
466 : ! **************************************************************************************************
467 128 : SUBROUTINE vwn_lsd_eval(rho_set, deriv_set, order, vwn_params)
468 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
469 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
470 : INTEGER, INTENT(in) :: order
471 : TYPE(section_vals_type), POINTER :: vwn_params
472 :
473 : CHARACTER(len=*), PARAMETER :: routineN = 'vwn_lsd_eval'
474 :
475 : INTEGER :: handle, npoints, TYPE
476 : INTEGER, DIMENSION(2, 3) :: bo
477 : REAL(KIND=dp) :: epsilon_rho, sc
478 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
479 32 : POINTER :: dummy, e_0, e_a, e_aa, e_aaa, e_aab, &
480 32 : e_ab, e_abb, e_b, e_bb, e_bbb, rhoa, &
481 32 : rhob
482 : TYPE(xc_derivative_type), POINTER :: deriv
483 :
484 32 : CALL timeset(routineN, handle)
485 :
486 32 : CALL section_vals_val_get(vwn_params, "scale_c", r_val=sc)
487 :
488 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
489 32 : local_bounds=bo, rho_cutoff=epsilon_rho)
490 32 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
491 32 : CALL vwn_init(epsilon_rho, vwn_params)
492 :
493 32 : dummy => rhoa
494 :
495 32 : e_0 => dummy
496 32 : e_a => dummy
497 32 : e_b => dummy
498 32 : e_aa => dummy
499 32 : e_bb => dummy
500 32 : e_ab => dummy
501 32 : e_aaa => dummy
502 32 : e_bbb => dummy
503 32 : e_aab => dummy
504 32 : e_abb => dummy
505 :
506 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
507 32 : allocate_deriv=.TRUE.)
508 32 : CALL xc_derivative_get(deriv, deriv_data=e_0)
509 32 : IF (order >= 1) THEN
510 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
511 32 : allocate_deriv=.TRUE.)
512 32 : CALL xc_derivative_get(deriv, deriv_data=e_a)
513 :
514 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
515 32 : allocate_deriv=.TRUE.)
516 32 : CALL xc_derivative_get(deriv, deriv_data=e_b)
517 : END IF
518 32 : IF (order >= 2) THEN
519 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
520 0 : allocate_deriv=.TRUE.)
521 0 : CALL xc_derivative_get(deriv, deriv_data=e_aa)
522 :
523 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
524 0 : allocate_deriv=.TRUE.)
525 0 : CALL xc_derivative_get(deriv, deriv_data=e_bb)
526 :
527 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
528 0 : allocate_deriv=.TRUE.)
529 0 : CALL xc_derivative_get(deriv, deriv_data=e_ab)
530 : END IF
531 32 : IF (order >= 3) THEN
532 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
533 0 : allocate_deriv=.TRUE.)
534 0 : CALL xc_derivative_get(deriv, deriv_data=e_aaa)
535 :
536 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
537 0 : allocate_deriv=.TRUE.)
538 0 : CALL xc_derivative_get(deriv, deriv_data=e_bbb)
539 :
540 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
541 0 : allocate_deriv=.TRUE.)
542 0 : CALL xc_derivative_get(deriv, deriv_data=e_aab)
543 :
544 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
545 0 : allocate_deriv=.TRUE.)
546 0 : CALL xc_derivative_get(deriv, deriv_data=e_abb)
547 :
548 : END IF
549 32 : IF (order > 3) THEN
550 0 : CPABORT("derivatives bigger than 3 not implemented")
551 : END IF
552 :
553 32 : CALL section_vals_val_get(vwn_params, "functional_type", i_val=TYPE)
554 30 : SELECT CASE (TYPE)
555 : CASE (do_vwn3)
556 :
557 : !$OMP PARALLEL DEFAULT(NONE) &
558 : !$OMP SHARED(rhoa, rhob, e_0, e_a, e_b, e_aa, e_bb, e_ab) &
559 : !$OMP SHARED(e_aaa, e_bbb, e_aab, e_abb, order, npoints) &
560 30 : !$OMP SHARED(sc)
561 :
562 : CALL vwn3_lsd_calc( &
563 : rhoa=rhoa, rhob=rhob, e_0=e_0, &
564 : e_a=e_a, e_b=e_b, &
565 : e_aa=e_aa, e_bb=e_bb, e_ab=e_ab, &
566 : e_aaa=e_aaa, e_bbb=e_bbb, e_aab=e_aab, e_abb=e_abb, &
567 : order=order, npoints=npoints, &
568 : sc=sc)
569 :
570 : !$OMP END PARALLEL
571 :
572 : CASE (do_vwn5)
573 :
574 : !$OMP PARALLEL DEFAULT(NONE) &
575 : !$OMP SHARED(rhoa, rhob, e_0, e_a, e_b, e_aa, e_bb, e_ab) &
576 : !$OMP SHARED(e_aaa, e_bbb, e_aab, e_abb, order, npoints) &
577 2 : !$OMP SHARED(sc)
578 :
579 : CALL vwn5_lsd_calc( &
580 : rhoa=rhoa, rhob=rhob, e_0=e_0, &
581 : e_a=e_a, e_b=e_b, &
582 : e_aa=e_aa, e_bb=e_bb, e_ab=e_ab, &
583 : e_aaa=e_aaa, e_bbb=e_bbb, e_aab=e_aab, e_abb=e_abb, &
584 : order=order, npoints=npoints, &
585 : sc=sc)
586 :
587 : !$OMP END PARALLEL
588 :
589 : CASE DEFAULT
590 32 : CPABORT(" Only functionals VWN3 and VWN5 are supported")
591 : END SELECT
592 :
593 32 : CALL timestop(handle)
594 :
595 32 : END SUBROUTINE vwn_lsd_eval
596 :
597 : ! **************************************************************************************************
598 : !> \brief ...
599 : !> \param rhoa ...
600 : !> \param rhob ...
601 : !> \param e_0 ...
602 : !> \param e_a ...
603 : !> \param e_b ...
604 : !> \param e_aa ...
605 : !> \param e_bb ...
606 : !> \param e_ab ...
607 : !> \param e_aaa ...
608 : !> \param e_bbb ...
609 : !> \param e_aab ...
610 : !> \param e_abb ...
611 : !> \param order ...
612 : !> \param npoints ...
613 : !> \param sc ...
614 : ! **************************************************************************************************
615 30 : SUBROUTINE vwn3_lsd_calc(rhoa, rhob, e_0, e_a, e_b, e_aa, e_bb, e_ab, &
616 : e_aaa, e_bbb, e_aab, e_abb, order, npoints, sc)
617 :
618 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: rhoa, rhob, e_0, e_a, e_b, e_aa, e_bb, &
619 : e_ab, e_aaa, e_bbb, e_aab, e_abb
620 : INTEGER, INTENT(in) :: order, npoints
621 : REAL(KIND=dp) :: sc
622 :
623 : INTEGER :: ip
624 : REAL(KIND=dp) :: ap, bp, cp, myrho, myrhoa, myrhob, Qf, Qp, t1, t10, t100, t101, t1019, &
625 : t102, t1031, t105, t108, t109, t11, t110, t113, t114, t115, t116, t117, t119, t120, t121, &
626 : t124, t125, t126, t127, t130, t132, t133, t134, t136, t137, t14, t144, t145, t146, t147, &
627 : t15, t150, t151, t154, t155, t156, t159, t16, t160, t161, t162, t165, t168, t169, t172, &
628 : t173, t174, t175, t176, t178, t179, t180, t183, t184, t187, t189, t19, t190, t191, t193, &
629 : t194, t2, t202, t203, t206, t209, t21, t212, t213, t216, t217, t218, t219, t220, t221, &
630 : t222, t223, t225, t226, t230, t231, t235, t236, t237, t24, t240
631 : REAL(KIND=dp) :: t241, t242, t244, t245, t246, t25, t251, t254, t255, t258, t259, t26, t260, &
632 : t261, t267, t268, t269, t27, t270, t273, t276, t279, t281, t282, t283, t284, t286, t287, &
633 : t29, t291, t292, t295, t296, t299, t3, t302, t306, t307, t309, t310, t311, t315, t316, &
634 : t319, t32, t324, t325, t332, t333, t334, t337, t339, t342, t343, t346, t349, t350, t351, &
635 : t352, t353, t354, t356, t357, t36, t363, t364, t365, t368, t373, t376, t377, t38, t380, &
636 : t381, t382, t388, t389, t390, t391, t394, t397, t4, t400, t402, t403, t404, t405, t407, &
637 : t408, t412, t413, t42, t420, t424, t425, t427, t428, t429, t43
638 : REAL(KIND=dp) :: t433, t434, t437, t44, t442, t443, t45, t451, t452, t455, t457, t458, t46, &
639 : t461, t464, t467, t470, t471, t472, t475, t479, t48, t482, t485, t488, t489, t49, t490, &
640 : t494, t497, t499, t5, t500, t501, t504, t505, t508, t509, t51, t510, t513, t516, t519, &
641 : t522, t525, t526, t528, t531, t536, t538, t54, t540, t541, t549, t55, t559, t563, t565, &
642 : t567, t570, t573, t58, t583, t589, t59, t595, t596, t6, t601, t603, t613, t62, t625, t64, &
643 : t665, t67, t68, t69, t690, t693, t696, t7, t705, t71, t710, t719, t720, t721, t727, t729, &
644 : t732, t74, t743, t745, t748, t752, t755, t758, t769, t78
645 : REAL(KIND=dp) :: t792, t793, t798, t80, t800, t804, t807, t808, t810, t814, t820, t823, &
646 : t835, t836, t838, t841, t85, t851, t853, t86, t860, t863, t872, t879, t88, t89, t90, t91, &
647 : t92, t947, t95, t950, t953, t956, t96, t967, t97, t98, t984, t986, t990, x0p
648 :
649 30 : cp = c
650 30 : ap = a
651 30 : bp = b
652 30 : x0p = x0
653 30 : Qp = SQRT(4.0_dp*cp - bp*bp)
654 30 : Qf = SQRT(4.0_dp*cf - bf*bf)
655 :
656 30 : !$OMP DO
657 : DO ip = 1, npoints
658 708000 : myrhoa = MAX(rhoa(ip), 0.0_dp)
659 708000 : myrhob = MAX(rhob(ip), 0.0_dp)
660 708000 : myrho = myrhoa + myrhob
661 708000 : IF (myrho > eps_rho) THEN
662 694537 : myrhoa = MAX(EPSILON(0.0_dp)*1.e4_dp, myrhoa)
663 694537 : myrhob = MAX(EPSILON(0.0_dp)*1.e4_dp, myrhob)
664 694537 : myrho = myrhoa + myrhob
665 :
666 694537 : t1 = 0.1e1_dp/0.3141592654e1_dp
667 694537 : t2 = myrho
668 694537 : t3 = 0.1e1_dp/t2
669 694537 : t4 = t1*t3
670 694537 : t5 = t4**0.3333333332e0_dp
671 694537 : t6 = 0.9085602964e0_dp*t5
672 694537 : t7 = t4**0.1666666666e0_dp
673 694537 : t10 = t6 + 0.9531842930e0_dp*bp*t7 + cp
674 694537 : t11 = 0.1e1_dp/t10
675 694537 : t14 = LOG(0.9085602964e0_dp*t5*t11)
676 694537 : t15 = 0.1906368586e1_dp*t7
677 694537 : t16 = t15 + bp
678 694537 : t19 = ATAN(Qp/t16)
679 694537 : t21 = 0.1e1_dp/Qp
680 694537 : t24 = bp*x0p
681 694537 : t25 = 0.9531842930e0_dp*t7
682 694537 : t26 = t25 - x0p
683 694537 : t27 = t26**0.20e1_dp
684 694537 : t29 = LOG(t27*t11)
685 694537 : t32 = 0.20e1_dp*bp + 0.40e1_dp*x0p
686 694537 : t36 = x0p**0.20e1_dp
687 694537 : t38 = 0.1e1_dp/(t36 + t24 + cp)
688 : t42 = ap*(t14 + 0.20e1_dp*bp*t19*t21 - t24*(t29 + t32*t19 &
689 694537 : *t21)*t38)
690 694537 : t43 = myrhoa - myrhob
691 694537 : t44 = t43*t3
692 694537 : t45 = 0.10e1_dp + t44
693 694537 : t46 = t45**0.1333333333e1_dp
694 694537 : t48 = 0.10e1_dp - t44
695 694537 : t49 = t48**0.1333333333e1_dp
696 694537 : t51 = 0.1923661050e1_dp*t46 + 0.1923661050e1_dp*t49 - 0.3847322100e1_dp
697 694537 : t54 = t6 + 0.9531842930e0_dp*bf*t7 + cf
698 694537 : t55 = 0.1e1_dp/t54
699 694537 : t58 = LOG(0.9085602964e0_dp*t5*t55)
700 694537 : t59 = t15 + bf
701 694537 : t62 = ATAN(Qf/t59)
702 694537 : t64 = 0.1e1_dp/Qf
703 694537 : t67 = bf*x0f
704 694537 : t68 = t25 - x0f
705 694537 : t69 = t68**0.20e1_dp
706 694537 : t71 = LOG(t69*t55)
707 694537 : t74 = 0.20e1_dp*bf + 0.40e1_dp*x0f
708 694537 : t78 = x0f**0.20e1_dp
709 694537 : t80 = 0.1e1_dp/(t78 + t67 + cf)
710 : t85 = af*(t58 + 0.20e1_dp*bf*t62*t64 - t67*(t71 + t74*t62 &
711 694537 : *t64)*t80) - t42
712 694537 : t86 = t51*t85
713 :
714 694537 : e_0(ip) = e_0(ip) + ((t42 + t86)*t2)*sc
715 :
716 694537 : IF (order >= 1) THEN
717 694537 : t88 = t4**(-0.6666666668e0_dp)
718 694537 : t89 = t88*t11
719 694537 : t90 = t2**2
720 694537 : t91 = 0.1e1_dp/t90
721 694537 : t92 = t1*t91
722 694537 : t95 = t10**2
723 694537 : t96 = 0.1e1_dp/t95
724 694537 : t97 = t5*t96
725 694537 : t98 = t88*t1
726 694537 : t100 = 0.3028534320e0_dp*t98*t91
727 694537 : t101 = t4**(-0.8333333334e0_dp)
728 694537 : t102 = bp*t101
729 694537 : t105 = -t100 - 0.1588640488e0_dp*t102*t92
730 694537 : t108 = -0.3028534320e0_dp*t89*t92 - 0.9085602964e0_dp*t97*t105
731 694537 : t109 = t4**(-0.3333333332e0_dp)
732 694537 : t110 = t108*t109
733 694537 : t113 = t16**2
734 694537 : t114 = 0.1e1_dp/t113
735 694537 : t115 = bp*t114
736 694537 : t116 = t115*t101
737 694537 : t117 = Qp**2
738 694537 : t119 = 0.1e1_dp + t117*t114
739 694537 : t120 = 0.1e1_dp/t119
740 694537 : t121 = t92*t120
741 694537 : t124 = t26**0.10e1_dp
742 694537 : t125 = t124*t11
743 694537 : t126 = t101*t1
744 694537 : t127 = t126*t91
745 694537 : t130 = t27*t96
746 694537 : t132 = -0.3177280976e0_dp*t125*t127 - t130*t105
747 694537 : t133 = t26**(-0.20e1_dp)
748 694537 : t134 = t132*t133
749 694537 : t136 = t32*t114
750 694537 : t137 = t136*t101
751 : t144 = ap*(0.1100642416e1_dp*t110*t10 + 0.6354561950e0_dp*t116* &
752 694537 : t121 - t24*(t134*t10 + 0.3177280975e0_dp*t137*t121)*t38)
753 694537 : t145 = t45**0.333333333e0_dp
754 694537 : t146 = t43*t91
755 694537 : t147 = t3 - t146
756 694537 : t150 = t48**0.333333333e0_dp
757 694537 : t151 = -t147
758 694537 : t154 = 0.2564881399e1_dp*t145*t147 + 0.2564881399e1_dp*t150*t151
759 694537 : t155 = t154*t85
760 694537 : t156 = t88*t55
761 694537 : t159 = t54**2
762 694537 : t160 = 0.1e1_dp/t159
763 694537 : t161 = t5*t160
764 694537 : t162 = bf*t101
765 694537 : t165 = -t100 - 0.1588640488e0_dp*t162*t92
766 694537 : t168 = -0.3028534320e0_dp*t156*t92 - 0.9085602964e0_dp*t161*t165
767 694537 : t169 = t168*t109
768 694537 : t172 = t59**2
769 694537 : t173 = 0.1e1_dp/t172
770 694537 : t174 = bf*t173
771 694537 : t175 = t174*t101
772 694537 : t176 = Qf**2
773 694537 : t178 = 0.1e1_dp + t176*t173
774 694537 : t179 = 0.1e1_dp/t178
775 694537 : t180 = t92*t179
776 694537 : t183 = t68**0.10e1_dp
777 694537 : t184 = t183*t55
778 694537 : t187 = t69*t160
779 694537 : t189 = -0.3177280976e0_dp*t184*t127 - t187*t165
780 694537 : t190 = t68**(-0.20e1_dp)
781 694537 : t191 = t189*t190
782 694537 : t193 = t74*t173
783 694537 : t194 = t193*t101
784 : t202 = af*(0.1100642416e1_dp*t169*t54 + 0.6354561950e0_dp*t175* &
785 : t180 - t67*(t191*t54 + 0.3177280975e0_dp*t194*t180)*t80) - &
786 694537 : t144
787 694537 : t203 = t51*t202
788 :
789 694537 : e_a(ip) = e_a(ip) + ((t144 + t155 + t203)*t2 + t42 + t86)*sc
790 :
791 694537 : t206 = -t3 - t146
792 694537 : t209 = -t206
793 694537 : t212 = 0.2564881399e1_dp*t145*t206 + 0.2564881399e1_dp*t150*t209
794 694537 : t213 = t212*t85
795 :
796 694537 : e_b(ip) = e_b(ip) + ((t144 + t213 + t203)*t2 + t42 + t86)*sc
797 :
798 : END IF
799 694537 : IF (order >= 2) THEN
800 0 : t216 = t4**(-0.1666666667e1_dp)
801 0 : t217 = t216*t11
802 0 : t218 = 0.3141592654e1_dp**2
803 0 : t219 = 0.1e1_dp/t218
804 0 : t220 = t90**2
805 0 : t221 = 0.1e1_dp/t220
806 0 : t222 = t219*t221
807 0 : t223 = t217*t222
808 0 : t225 = t88*t96
809 0 : t226 = t92*t105
810 0 : t230 = 0.1e1_dp/t90/t2
811 0 : t231 = t1*t230
812 0 : t235 = 0.1e1_dp/t95/t10
813 0 : t236 = t5*t235
814 0 : t237 = t105**2
815 0 : t240 = t216*t219
816 0 : t241 = t240*t221
817 0 : t242 = 0.2019022880e0_dp*t241
818 0 : t244 = 0.6057068640e0_dp*t98*t230
819 0 : t245 = t4**(-0.1833333333e1_dp)
820 0 : t246 = bp*t245
821 : t251 = -t242 + t244 - 0.1323867073e0_dp*t246*t222 + 0.3177280976e0_dp &
822 0 : *t102*t231
823 : t254 = -0.2019022880e0_dp*t223 + 0.6057068640e0_dp*t225*t226 + 0.6057068640e0_dp &
824 : *t89*t231 + 0.1817120593e1_dp*t236*t237 - 0.9085602964e0_dp &
825 0 : *t97*t251
826 0 : t255 = t254*t109
827 0 : t258 = t4**(-0.1333333333e1_dp)
828 0 : t259 = t108*t258
829 0 : t260 = t10*t1
830 0 : t261 = t260*t91
831 0 : t267 = 0.1e1_dp/t113/t16
832 0 : t268 = bp*t267
833 0 : t269 = t268*t216
834 0 : t270 = t222*t120
835 0 : t273 = t115*t245
836 0 : t276 = t231*t120
837 0 : t279 = t113**2
838 0 : t281 = 0.1e1_dp/t279/t16
839 0 : t282 = bp*t281
840 0 : t283 = t282*t216
841 0 : t284 = t119**2
842 0 : t286 = 0.1e1_dp/t284*t117
843 0 : t287 = t222*t286
844 0 : t291 = t124*t96
845 0 : t292 = t291*t101
846 0 : t295 = t245*t219
847 0 : t296 = t295*t221
848 0 : t299 = t126*t230
849 0 : t302 = t27*t235
850 : t306 = 0.5047557200e-1_dp*t223 + 0.6354561952e0_dp*t292*t226 - 0.2647734147e0_dp &
851 : *t125*t296 + 0.6354561952e0_dp*t125*t299 + 0.2e1_dp* &
852 0 : t302*t237 - t130*t251
853 0 : t307 = t306*t133
854 0 : t309 = t26**(-0.30e1_dp)
855 0 : t310 = t132*t309
856 0 : t311 = t310*t10
857 0 : t315 = t32*t267
858 0 : t316 = t315*t216
859 0 : t319 = t136*t245
860 0 : t324 = t32*t281
861 0 : t325 = t324*t216
862 : t332 = ap*(0.1100642416e1_dp*t255*t10 + 0.3668808052e0_dp*t259* &
863 : t261 + 0.1100642416e1_dp*t110*t105 + 0.4038045758e0_dp*t269*t270 &
864 : + 0.5295468292e0_dp*t273*t270 - 0.1270912390e1_dp*t116*t276 - 0.4038045758e0_dp &
865 : *t283*t287 - t24*(t307*t10 + 0.3177280976e0_dp* &
866 : t311*t127 + t134*t105 + 0.2019022879e0_dp*t316*t270 + 0.2647734146e0_dp &
867 : *t319*t270 - 0.6354561950e0_dp*t137*t276 - 0.2019022879e0_dp &
868 0 : *t325*t287)*t38)
869 0 : t333 = t45**(-0.666666667e0_dp)
870 0 : t334 = t147**2
871 0 : t337 = t43*t230
872 0 : t339 = -0.2e1_dp*t91 + 0.2e1_dp*t337
873 0 : t342 = t48**(-0.666666667e0_dp)
874 0 : t343 = t151**2
875 0 : t346 = -t339
876 : t349 = 0.8549604655e0_dp*t333*t334 + 0.2564881399e1_dp*t145*t339 &
877 0 : + 0.8549604655e0_dp*t342*t343 + 0.2564881399e1_dp*t150*t346
878 0 : t350 = t349*t85
879 0 : t351 = t154*t202
880 0 : t352 = 0.2e1_dp*t351
881 0 : t353 = t216*t55
882 0 : t354 = t353*t222
883 0 : t356 = t88*t160
884 0 : t357 = t92*t165
885 0 : t363 = 0.1e1_dp/t159/t54
886 0 : t364 = t5*t363
887 0 : t365 = t165**2
888 0 : t368 = bf*t245
889 : t373 = -t242 + t244 - 0.1323867073e0_dp*t368*t222 + 0.3177280976e0_dp &
890 0 : *t162*t231
891 : t376 = -0.2019022880e0_dp*t354 + 0.6057068640e0_dp*t356*t357 + 0.6057068640e0_dp &
892 : *t156*t231 + 0.1817120593e1_dp*t364*t365 - 0.9085602964e0_dp &
893 0 : *t161*t373
894 0 : t377 = t376*t109
895 0 : t380 = t168*t258
896 0 : t381 = t54*t1
897 0 : t382 = t381*t91
898 0 : t388 = 0.1e1_dp/t172/t59
899 0 : t389 = bf*t388
900 0 : t390 = t389*t216
901 0 : t391 = t222*t179
902 0 : t394 = t174*t245
903 0 : t397 = t231*t179
904 0 : t400 = t172**2
905 0 : t402 = 0.1e1_dp/t400/t59
906 0 : t403 = bf*t402
907 0 : t404 = t403*t216
908 0 : t405 = t178**2
909 0 : t407 = 0.1e1_dp/t405*t176
910 0 : t408 = t222*t407
911 0 : t412 = t183*t160
912 0 : t413 = t412*t101
913 0 : t420 = t69*t363
914 : t424 = 0.5047557200e-1_dp*t354 + 0.6354561952e0_dp*t413*t357 - 0.2647734147e0_dp &
915 : *t184*t296 + 0.6354561952e0_dp*t184*t299 + 0.2e1_dp* &
916 0 : t420*t365 - t187*t373
917 0 : t425 = t424*t190
918 0 : t427 = t68**(-0.30e1_dp)
919 0 : t428 = t189*t427
920 0 : t429 = t428*t54
921 0 : t433 = t74*t388
922 0 : t434 = t433*t216
923 0 : t437 = t193*t245
924 0 : t442 = t74*t402
925 0 : t443 = t442*t216
926 : t451 = af*(0.1100642416e1_dp*t377*t54 + 0.3668808052e0_dp*t380* &
927 : t382 + 0.1100642416e1_dp*t169*t165 + 0.4038045758e0_dp*t390*t391 &
928 : + 0.5295468292e0_dp*t394*t391 - 0.1270912390e1_dp*t175*t397 - 0.4038045758e0_dp &
929 : *t404*t408 - t67*(t425*t54 + 0.3177280976e0_dp* &
930 : t429*t127 + t191*t165 + 0.2019022879e0_dp*t434*t391 + 0.2647734146e0_dp &
931 : *t437*t391 - 0.6354561950e0_dp*t194*t397 - 0.2019022879e0_dp &
932 0 : *t443*t408)*t80) - t332
933 0 : t452 = t51*t451
934 0 : t455 = 0.2e1_dp*t144
935 0 : t457 = 0.2e1_dp*t203
936 :
937 0 : e_aa(ip) = e_aa(ip) + ((t332 + t350 + t352 + t452)*t2 + t455 + 0.2e1_dp*t155 + t457)*sc
938 :
939 0 : t458 = t333*t147
940 0 : t461 = t145*t43
941 0 : t464 = t342*t151
942 0 : t467 = t150*t43
943 : t470 = 0.8549604655e0_dp*t458*t206 + 0.5129762798e1_dp*t461*t230 &
944 0 : + 0.8549604655e0_dp*t464*t209 - 0.5129762798e1_dp*t467*t230
945 0 : t471 = t470*t85
946 0 : t472 = t212*t202
947 :
948 : e_ab(ip) = e_ab(ip) + ((t332 + t471 + t351 + t472 + t452)*t2 + t455 + t155 + t457 &
949 0 : + t213)*sc
950 0 : t475 = t206**2
951 0 : t479 = 0.2e1_dp*t91 + 0.2e1_dp*t337
952 0 : t482 = t209**2
953 0 : t485 = -t479
954 : t488 = 0.8549604655e0_dp*t333*t475 + 0.2564881399e1_dp*t145*t479 &
955 0 : + 0.8549604655e0_dp*t342*t482 + 0.2564881399e1_dp*t150*t485
956 0 : t489 = t488*t85
957 0 : t490 = 0.2e1_dp*t472
958 :
959 0 : e_bb(ip) = e_bb(ip) + ((t332 + t489 + t490 + t452)*t2 + t455 + 0.2e1_dp*t213 + t457)*sc
960 :
961 : END IF
962 694537 : IF (order >= 3) THEN
963 0 : t494 = t4**(-0.2666666667e1_dp)
964 0 : t497 = 0.1e1_dp/t218/0.3141592654e1_dp
965 0 : t499 = 0.1e1_dp/t220/t90
966 0 : t500 = t497*t499
967 0 : t501 = t494*t11*t500
968 0 : t504 = t222*t105
969 0 : t505 = t216*t96*t504
970 0 : t508 = 0.1e1_dp/t220/t2
971 0 : t509 = t219*t508
972 0 : t510 = t217*t509
973 0 : t513 = t92*t237
974 0 : t516 = t231*t105
975 0 : t519 = t92*t251
976 0 : t522 = t1*t221
977 0 : t525 = t95**2
978 0 : t526 = 0.1e1_dp/t525
979 0 : t528 = t237*t105
980 0 : t531 = t105*t251
981 0 : t536 = 0.3365038134e0_dp*t494*t497*t499
982 0 : t538 = 0.1211413728e1_dp*t240*t508
983 0 : t540 = 0.1817120592e1_dp*t98*t221
984 0 : t541 = t4**(-0.2833333333e1_dp)
985 : t549 = -t536 + t538 - t540 - 0.2427089633e0_dp*bp*t541*t500 + 0.7943202439e0_dp &
986 0 : *t246*t509 - 0.9531842928e0_dp*t102*t522
987 0 : t559 = t500*t120
988 0 : t563 = 0.1e1_dp/t279/t113
989 0 : t565 = t4**(-0.2500000000e1_dp)
990 0 : t567 = t500*t286
991 0 : t570 = t509*t120
992 0 : t573 = t4**(-0.2666666666e1_dp)
993 0 : t583 = t522*t120
994 0 : t589 = t509*t286
995 0 : t595 = t279**2
996 0 : t596 = 0.1e1_dp/t595
997 0 : t601 = t117**2
998 0 : t603 = t500/t284/t119*t601
999 0 : t613 = t4**(-0.2333333333e1_dp)
1000 0 : t625 = 0.1e1_dp/t279
1001 0 : t665 = t26**(-0.40e1_dp)
1002 0 : t690 = t541*t497*t499
1003 0 : t693 = t295*t508
1004 0 : t696 = t126*t221
1005 : t705 = 0.8412595335e-1_dp*t501 - 0.1514267160e0_dp*t505 - 0.3028534320e0_dp &
1006 : *t510 - 0.1906368585e1_dp*t124*t235*t101*t513 + 0.7943202441e0_dp &
1007 : *t291*t245*t504 - 0.1906368585e1_dp*t292*t516 + 0.9531842928e0_dp &
1008 : *t292*t519 + 0.4206297667e-1_dp*t11*t573*t500 - 0.4854179269e0_dp &
1009 : *t125*t690 + 0.1588640488e1_dp*t125*t693 - 0.1906368586e1_dp &
1010 : *t125*t696 - 0.6e1_dp*t27*t526*t528 + 0.6e1_dp*t302 &
1011 0 : *t531 - t130*t549
1012 : t710 = 0.6354561952e0_dp*t306*t309*t10*t127 - 0.1211413727e1_dp* &
1013 : t316*t570 + 0.1924500894e0_dp*t32*t625*t565*t559 + 0.3365038132e0_dp &
1014 : *t315*t494*t559 - 0.4490502088e0_dp*t32*t563*t565 &
1015 : *t567 - 0.1588640487e1_dp*t319*t570 + 0.1682519066e0_dp*t315*t573 &
1016 : *t559 + 0.4854179267e0_dp*t136*t541*t559 - 0.1682519066e0_dp* &
1017 : t324*t573*t567 + 0.1906368585e1_dp*t137*t583 + 0.1211413727e1_dp &
1018 : *t325*t589 - 0.3365038132e0_dp*t324*t494*t567 + 0.2566001193e0_dp &
1019 : *t32*t596*t565*t603 + t134*t251 + 0.6354561952e0_dp*t310 &
1020 : *t105*t127 - 0.6354561952e0_dp*t311*t299 + 0.1514267160e0_dp &
1021 : *t132*t665*t10*t241 + 0.2647734147e0_dp*t311*t296 + t705* &
1022 0 : t133*t10 + 0.2e1_dp*t307*t105
1023 : t719 = 0.1100642416e1_dp*(-0.3365038134e0_dp*t501 + 0.6057068641e0_dp* &
1024 : t505 + 0.1211413728e1_dp*t510 - 0.1817120592e1_dp*t88*t235*t513 &
1025 : - 0.1817120592e1_dp*t225*t516 + 0.9085602960e0_dp*t225*t519 - 0.1817120592e1_dp &
1026 : *t89*t522 - 0.5451361779e1_dp*t5*t526*t528 + 0.5451361779e1_dp &
1027 : *t236*t531 - 0.9085602964e0_dp*t97*t549)*t109* &
1028 : t10 + 0.2201284832e1_dp*t255*t105 + 0.6730076265e0_dp*t268*t494 &
1029 : *t559 - 0.8981004177e0_dp*bp*t563*t565*t567 - 0.3177280975e1_dp &
1030 : *t273*t570 + 0.3365038132e0_dp*t268*t573*t559 + 0.9708358534e0_dp &
1031 : *t115*t541*t559 - 0.3365038132e0_dp*t282*t573*t567 + &
1032 : 0.3812737170e1_dp*t116*t583 + 0.7337616104e0_dp*t254*t258*t261 &
1033 : + 0.2422827454e1_dp*t283*t589 - 0.6730076265e0_dp*t282*t494*t567 &
1034 : + 0.5132002385e0_dp*bp*t596*t565*t603 + 0.7337616104e0_dp* &
1035 : t259*t226 + 0.1100642416e1_dp*t110*t251 - 0.7337616104e0_dp*t259 &
1036 : *t260*t230 + 0.4891744068e0_dp*t108*t613*t10*t219*t221 &
1037 : - t24*t710*t38 - 0.2422827454e1_dp*t269*t570 + 0.3849001789e0_dp &
1038 0 : *bp*t625*t565*t559
1039 0 : t720 = ap*t719
1040 0 : t721 = t45**(-0.1666666667e1_dp)
1041 0 : t727 = t43*t221
1042 0 : t729 = 0.6e1_dp*t230 - 0.6e1_dp*t727
1043 0 : t732 = t48**(-0.1666666667e1_dp)
1044 0 : t743 = t349*t202
1045 0 : t745 = t154*t451
1046 0 : t748 = t500*t179
1047 0 : t752 = t500*t407
1048 0 : t755 = t522*t179
1049 0 : t758 = t509*t407
1050 0 : t769 = t68**(-0.40e1_dp)
1051 0 : t792 = t400**2
1052 0 : t793 = 0.1e1_dp/t792
1053 0 : t798 = t176**2
1054 0 : t800 = t500/t405/t178*t798
1055 0 : t804 = t494*t55*t500
1056 0 : t807 = t222*t165
1057 0 : t808 = t216*t160*t807
1058 0 : t810 = t353*t509
1059 0 : t814 = t92*t365
1060 0 : t820 = t231*t165
1061 0 : t823 = t92*t373
1062 0 : t835 = t159**2
1063 0 : t836 = 0.1e1_dp/t835
1064 0 : t838 = t365*t165
1065 0 : t841 = t165*t373
1066 : t851 = -t536 + t538 - t540 - 0.2427089633e0_dp*bf*t541*t500 + 0.7943202439e0_dp &
1067 0 : *t368*t509 - 0.9531842928e0_dp*t162*t522
1068 : t853 = 0.8412595335e-1_dp*t804 - 0.1514267160e0_dp*t808 - 0.3028534320e0_dp &
1069 : *t810 - 0.1906368585e1_dp*t183*t363*t101*t814 + 0.7943202441e0_dp &
1070 : *t412*t245*t807 - 0.1906368585e1_dp*t413*t820 + 0.9531842928e0_dp &
1071 : *t413*t823 + 0.4206297667e-1_dp*t55*t573*t500 - 0.4854179269e0_dp &
1072 : *t184*t690 + 0.1588640488e1_dp*t184*t693 - 0.1906368586e1_dp &
1073 : *t184*t696 - 0.6e1_dp*t69*t836*t838 + 0.6e1_dp*t420 &
1074 0 : *t841 - t187*t851
1075 0 : t860 = t509*t179
1076 0 : t863 = 0.1e1_dp/t400
1077 0 : t872 = 0.1e1_dp/t400/t172
1078 : t879 = 0.2e1_dp*t425*t165 + t191*t373 + 0.6354561952e0_dp*t428* &
1079 : t165*t127 - 0.6354561952e0_dp*t429*t299 + 0.1514267160e0_dp*t189 &
1080 : *t769*t54*t241 + 0.2647734147e0_dp*t429*t296 + 0.1682519066e0_dp &
1081 : *t433*t573*t748 + 0.4854179267e0_dp*t193*t541*t748 - 0.1682519066e0_dp &
1082 : *t442*t573*t752 + 0.1906368585e1_dp*t194*t755 + &
1083 : 0.1211413727e1_dp*t443*t758 - 0.3365038132e0_dp*t442*t494*t752 &
1084 : + 0.2566001193e0_dp*t74*t793*t565*t800 + t853*t190*t54 &
1085 : + 0.6354561952e0_dp*t424*t427*t54*t127 - 0.1211413727e1_dp*t434 &
1086 : *t860 + 0.1924500894e0_dp*t74*t863*t565*t748 + 0.3365038132e0_dp &
1087 : *t433*t494*t748 - 0.4490502088e0_dp*t74*t872*t565*t752 &
1088 0 : - 0.1588640487e1_dp*t437*t860
1089 : t947 = 0.9708358534e0_dp*t174*t541*t748 - 0.3365038132e0_dp*t403 &
1090 : *t573*t752 + 0.3812737170e1_dp*t175*t755 + 0.2422827454e1_dp*t404 &
1091 : *t758 - t67*t879*t80 - 0.6730076265e0_dp*t403*t494*t752 &
1092 : + 0.5132002385e0_dp*bf*t793*t565*t800 + 0.2201284832e1_dp*t377 &
1093 : *t165 + 0.1100642416e1_dp*(-0.3365038134e0_dp*t804 + 0.6057068641e0_dp &
1094 : *t808 + 0.1211413728e1_dp*t810 - 0.1817120592e1_dp*t88*t363* &
1095 : t814 - 0.1817120592e1_dp*t356*t820 + 0.9085602960e0_dp*t356*t823 &
1096 : - 0.1817120592e1_dp*t156*t522 - 0.5451361779e1_dp*t5*t836*t838 &
1097 : + 0.5451361779e1_dp*t364*t841 - 0.9085602964e0_dp*t161*t851)* &
1098 : t109*t54 + 0.7337616104e0_dp*t376*t258*t382 + 0.1100642416e1_dp &
1099 : *t169*t373 + 0.7337616104e0_dp*t380*t357 - 0.7337616104e0_dp*t380 &
1100 : *t381*t230 + 0.4891744068e0_dp*t168*t613*t54*t219*t221 &
1101 : - 0.2422827454e1_dp*t390*t860 + 0.3849001789e0_dp*bf*t863*t565 &
1102 : *t748 + 0.6730076265e0_dp*t389*t494*t748 - 0.8981004177e0_dp &
1103 : *bf*t872*t565*t752 - 0.3177280975e1_dp*t394*t860 + 0.3365038132e0_dp &
1104 0 : *t389*t573*t748
1105 0 : t950 = t51*(af*t947 - t720)
1106 0 : t953 = 0.3e1_dp*t332
1107 0 : t956 = 0.3e1_dp*t452
1108 :
1109 : e_aaa(ip) = e_aaa(ip) + ((t720 + (-0.5699736440e0_dp*t721*t334*t147 + 0.2564881396e1_dp &
1110 : *t458*t339 + 0.2564881399e1_dp*t145*t729 - 0.5699736440e0_dp* &
1111 : t732*t343*t151 + 0.2564881396e1_dp*t464*t346 - 0.2564881399e1_dp &
1112 : *t150*t729)*t85 + 0.3e1_dp*t743 + 0.3e1_dp*t745 + t950)*t2 &
1113 0 : + t953 + 0.3e1_dp*t350 + 0.6e1_dp*t351 + t956)*sc
1114 :
1115 0 : t967 = 0.2e1_dp*t230 - 0.6e1_dp*t727
1116 0 : t984 = 0.2e1_dp*t470*t202
1117 0 : t986 = t212*t451
1118 0 : t990 = 0.2e1_dp*t471
1119 :
1120 : e_aab(ip) = e_aab(ip) + ((t720 + (-0.5699736440e0_dp*t721*t334*t206 + 0.3419841862e1_dp &
1121 : *t458*t337 + 0.8549604655e0_dp*t333*t339*t206 + 0.2564881399e1_dp &
1122 : *t145*t967 - 0.5699736440e0_dp*t732*t343*t209 - 0.3419841862e1_dp &
1123 : *t464*t337 + 0.8549604655e0_dp*t342*t346*t209 - 0.2564881399e1_dp &
1124 : *t150*t967)*t85 + t743 + t984 + 0.2e1_dp*t745 + t986 &
1125 0 : + t950)*t2 + t953 + t350 + 0.4e1_dp*t351 + t956 + t990 + t490)*sc
1126 :
1127 0 : t1019 = t488*t202
1128 :
1129 : e_abb(ip) = e_abb(ip) + ((t720 + (-0.5699736440e0_dp*t721*t147*t475 + 0.3419841862e1_dp &
1130 : *t333*t43*t230*t206 + 0.8549604655e0_dp*t458*t479 - 0.5129762798e1_dp &
1131 : *t145*t230 - 0.1538928839e2_dp*t461*t221 - 0.5699736440e0_dp &
1132 : *t732*t151*t482 - 0.3419841862e1_dp*t342*t43*t230 &
1133 : *t209 + 0.8549604655e0_dp*t464*t485 + 0.5129762798e1_dp*t150*t230 &
1134 : + 0.1538928839e2_dp*t467*t221)*t85 + t984 + t745 + t1019 + 0.2e1_dp &
1135 : *t986 + t950)*t2 + t953 + t990 + t352 + 0.4e1_dp*t472 + t956 &
1136 0 : + t489)*sc
1137 :
1138 0 : t1031 = -0.6e1_dp*t230 - 0.6e1_dp*t727
1139 :
1140 : e_bbb(ip) = e_bbb(ip) + ((t720 + (-0.5699736440e0_dp*t721*t475*t206 + 0.2564881396e1_dp &
1141 : *t333*t206*t479 + 0.2564881399e1_dp*t145*t1031 - 0.5699736440e0_dp &
1142 : *t732*t482*t209 + 0.2564881396e1_dp*t342*t209*t485 &
1143 : - 0.2564881399e1_dp*t150*t1031)*t85 + 0.3e1_dp*t1019 + 0.3e1_dp*t986 &
1144 0 : + t950)*t2 + t953 + 0.3e1_dp*t489 + 0.6e1_dp*t472 + t956)*sc
1145 :
1146 : END IF
1147 : END IF
1148 : END DO
1149 :
1150 : !$OMP END DO
1151 :
1152 30 : END SUBROUTINE vwn3_lsd_calc
1153 :
1154 : ! **************************************************************************************************
1155 : !> \brief ...
1156 : !> \param rhoa ...
1157 : !> \param rhob ...
1158 : !> \param e_0 ...
1159 : !> \param e_a ...
1160 : !> \param e_b ...
1161 : !> \param e_aa ...
1162 : !> \param e_bb ...
1163 : !> \param e_ab ...
1164 : !> \param e_aaa ...
1165 : !> \param e_bbb ...
1166 : !> \param e_aab ...
1167 : !> \param e_abb ...
1168 : !> \param order ...
1169 : !> \param npoints ...
1170 : !> \param sc ...
1171 : ! **************************************************************************************************
1172 2 : SUBROUTINE vwn5_lsd_calc(rhoa, rhob, e_0, e_a, e_b, e_aa, e_bb, e_ab, &
1173 : e_aaa, e_bbb, e_aab, e_abb, order, npoints, sc)
1174 :
1175 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: rhoa, rhob, e_0, e_a, e_b, e_aa, e_bb, &
1176 : e_ab, e_aaa, e_bbb, e_aab, e_abb
1177 : INTEGER, INTENT(in) :: order, npoints
1178 : REAL(KIND=dp) :: sc
1179 :
1180 : INTEGER :: ip
1181 : REAL(KIND=dp) :: ap, bp, cp, d2f0, myrho, myrhoa, myrhob, Qa, Qf, Qp, t1, t100, t101, t1010, &
1182 : t1013, t1019, t1020, t1025, t1027, t1030, t104, t1043, t106, t1084, t1085, t1086, t1088, &
1183 : t1089, t109, t1094, t1096, t1099, t11, t110, t1108, t1109, t111, t1110, t1111, t1113, &
1184 : t1114, t1116, t1118, t1119, t1120, t1125, t113, t1137, t1142, t1145, t1148, t1149, t1151, &
1185 : t1154, t1157, t116, t1160, t1165, t1166, t1168, t1171, t1181, t12, t120, t1205, t1208, &
1186 : t1211, t1218, t122, t1221, t1235, t1238, t1244, t1245, t125, t1250, t1252, t126, t127, &
1187 : t1284, t129, t1299, t13, t130, t131, t132, t133, t1341, t1344
1188 : REAL(KIND=dp) :: t1346, t1358, t136, t1360, t1361, t1363, t1365, t1368, t137, t1371, t1372, &
1189 : t1374, t1377, t1380, t1383, t1388, t1389, t1391, t1394, t140, t1404, t141, t142, t143, &
1190 : t144, t1455, t1469, t147, t1477, t148, t1480, t1483, t149, t1490, t1493, t15, t150, &
1191 : t1507, t151, t1510, t1516, t1517, t1522, t1524, t1527, t154, t155, t156, t1567, t157, &
1192 : t1570, t1576, t1578, t1580, t1582, t1588, t159, t1590, t1593, t1599, t16, t160, t1602, &
1193 : t1603, t1604, t1605, t1606, t1609, t161, t1610, t1611, t1613, t1615, t1620, t1622, t1626, &
1194 : t1639, t164, t1653, t1654, t1657, t1659, t1664, t1666, t1668, t167
1195 : REAL(KIND=dp) :: t1671, t1672, t1675, t168, t169, t1690, t1694, t1696, t1697, t1698, t17, &
1196 : t1701, t172, t1726, t173, t1730, t174, t175, t1750, t1754, t176, t1764, t1765, t1768, &
1197 : t1775, t1776, t1778, t178, t1782, t1783, t1785, t1786, t1789, t179, t1795, t1797, t18, &
1198 : t180, t1804, t1826, t1829, t183, t1832, t1836, t1838, t1839, t184, t185, t1850, t1853, &
1199 : t1858, t186, t1868, t1870, t1872, t1884, t1886, t189, t1891, t1893, t1894, t1897, t19, &
1200 : t1901, t1904, t191, t192, t193, t1939, t1945, t1949, t195, t196, t1975, t1979, t1986, &
1201 : t1998, t1999, t2, t20, t2010, t2011, t2014, t2019, t202, t2025, t203
1202 : REAL(KIND=dp) :: t204, t2048, t205, t206, t207, t208, t211, t212, t213, t214, t217, t220, &
1203 : t221, t224, t225, t226, t227, t228, t23, t230, t231, t232, t235, t236, t239, t24, t241, &
1204 : t242, t243, t245, t246, t252, t253, t254, t255, t256, t257, t258, t259, t260, t263, t264, &
1205 : t265, t266, t269, t27, t272, t273, t276, t277, t278, t279, t28, t280, t282, t283, t284, &
1206 : t287, t288, t29, t291, t293, t294, t295, t297, t298, t3, t304, t305, t306, t309, t312, &
1207 : t315, t316, t317, t32, t320, t321, t322, t323, t324, t325, t326, t327, t328, t329, t332, &
1208 : t333, t337, t338, t34, t340, t343, t344, t347, t350, t351, t352
1209 : REAL(KIND=dp) :: t353, t355, t356, t357, t359, t362, t363, t364, t365, t366, t367, t368, &
1210 : t369, t37, t370, t371, t372, t373, t375, t376, t379, t38, t383, t384, t385, t388, t389, &
1211 : t39, t390, t392, t393, t394, t399, t4, t40, t402, t403, t406, t407, t408, t409, t415, &
1212 : t416, t417, t418, t42, t421, t424, t427, t429, t430, t431, t432, t434, t435, t439, t440, &
1213 : t443, t444, t447, t45, t450, t454, t455, t457, t458, t459, t463, t464, t467, t472, t473, &
1214 : t479, t480, t481, t482, t483, t484, t485, t486, t487, t488, t489, t49, t490, t491, t492, &
1215 : t493, t494, t496, t497, t5, t503, t504, t505, t508, t51, t513
1216 : REAL(KIND=dp) :: t516, t517, t520, t521, t522, t528, t529, t530, t531, t534, t537, t54, &
1217 : t540, t542, t543, t544, t545, t547, t548, t55, t552, t553, t56, t560, t564, t565, t567, &
1218 : t568, t569, t57, t573, t574, t577, t582, t583, t591, t592, t593, t594, t595, t596, t597, &
1219 : t598, t599, t6, t60, t600, t601, t602, t603, t604, t605, t606, t607, t608, t61, t610, &
1220 : t611, t617, t618, t619, t622, t627, t630, t631, t634, t635, t636, t64, t642, t643, t644, &
1221 : t645, t648, t65, t651, t654, t656, t657, t658, t659, t661, t662, t666, t667, t674, t678, &
1222 : t679, t68, t681, t682, t683, t687, t688, t691, t696, t697, t70
1223 : REAL(KIND=dp) :: t704, t705, t706, t709, t712, t715, t716, t719, t722, t725, t728, t729, &
1224 : t73, t730, t732, t733, t735, t74, t741, t742, t743, t744, t745, t746, t748, t75, t750, &
1225 : t751, t752, t753, t755, t756, t757, t758, t759, t762, t763, t764, t766, t767, t769, t77, &
1226 : t771, t772, t773, t774, t777, t778, t779, t780, t782, t783, t784, t786, t789, t793, t796, &
1227 : t799, t8, t80, t802, t803, t804, t806, t808, t811, t812, t813, t814, t815, t816, t817, &
1228 : t818, t819, t820, t821, t822, t823, t824, t825, t826, t827, t828, t829, t830, t831, t832, &
1229 : t833, t834, t835, t84, t842, t851, t857, t859, t86, t862, t864
1230 : REAL(KIND=dp) :: t865, t866, t869, t870, t873, t874, t875, t878, t881, t884, t887, t89, &
1231 : t890, t891, t893, t896, t9, t90, t901, t903, t905, t906, t91, t914, t92, t93, t937, t942, &
1232 : t945, t948, t957, t96, t97, t972, t979, t982, t984, t986, t993, t996, x0p
1233 :
1234 2 : cp = c
1235 2 : ap = a
1236 2 : bp = b
1237 2 : x0p = x0
1238 2 : Qp = SQRT(4.0_dp*cp - bp*bp)
1239 2 : Qa = SQRT(4.0_dp*ca - ba*ba)
1240 2 : Qf = SQRT(4.0_dp*cf - bf*bf)
1241 2 : d2f0 = 4.0_dp/(9.0_dp*(2.0_dp**(1.0_dp/3.0_dp) - 1.0_dp))
1242 :
1243 2 : !$OMP DO
1244 :
1245 : DO ip = 1, npoints
1246 800 : myrhoa = MAX(rhoa(ip), 0.0_dp)
1247 800 : myrhob = MAX(rhob(ip), 0.0_dp)
1248 800 : myrho = myrhoa + myrhob
1249 800 : IF (myrho > eps_rho) THEN
1250 634 : myrhoa = MAX(EPSILON(0.0_dp)*1.e4_dp, myrhoa)
1251 634 : myrhob = MAX(EPSILON(0.0_dp)*1.e4_dp, myrhob)
1252 634 : myrho = myrhoa + myrhob
1253 :
1254 634 : t1 = myrhoa - myrhob
1255 634 : t2 = myrho
1256 634 : t3 = 0.1e1_dp/t2
1257 634 : t4 = t1*t3
1258 634 : t5 = 0.10e1_dp + t4
1259 634 : t6 = t5**0.1333333333e1_dp
1260 634 : t8 = 0.10e1_dp - t4
1261 634 : t9 = t8**0.1333333333e1_dp
1262 634 : t11 = 0.1923661050e1_dp*t6 + 0.1923661050e1_dp*t9 - 0.3847322100e1_dp
1263 634 : t12 = t4**0.40e1_dp
1264 634 : t13 = t11*t12
1265 634 : t15 = (0.10e1_dp - t13)*ap
1266 634 : t16 = 0.1e1_dp/0.3141592654e1_dp
1267 634 : t17 = t16*t3
1268 634 : t18 = t17**0.3333333332e0_dp
1269 634 : t19 = 0.9085602964e0_dp*t18
1270 634 : t20 = t17**0.1666666666e0_dp
1271 634 : t23 = t19 + 0.9531842930e0_dp*bp*t20 + cp
1272 634 : t24 = 0.1e1_dp/t23
1273 634 : t27 = LOG(0.9085602964e0_dp*t18*t24)
1274 634 : t28 = 0.1906368586e1_dp*t20
1275 634 : t29 = t28 + bp
1276 634 : t32 = ATAN(Qp/t29)
1277 634 : t34 = 0.1e1_dp/Qp
1278 634 : t37 = bp*x0p
1279 634 : t38 = 0.9531842930e0_dp*t20
1280 634 : t39 = t38 - x0p
1281 634 : t40 = t39**0.20e1_dp
1282 634 : t42 = LOG(t40*t24)
1283 634 : t45 = 0.20e1_dp*bp + 0.40e1_dp*x0p
1284 634 : t49 = x0p**0.20e1_dp
1285 634 : t51 = 0.1e1_dp/(t49 + t37 + cp)
1286 : t54 = t27 + 0.20e1_dp*bp*t32*t34 - t37*(t42 + t45*t32*t34) &
1287 634 : *t51
1288 634 : t55 = t15*t54
1289 634 : t56 = 0.10e1_dp - t12
1290 634 : t57 = t11*t56
1291 634 : t60 = t19 + 0.9531842930e0_dp*ba*t20 + ca
1292 634 : t61 = 0.1e1_dp/t60
1293 634 : t64 = LOG(0.9085602964e0_dp*t18*t61)
1294 634 : t65 = t28 + ba
1295 634 : t68 = ATAN(Qa/t65)
1296 634 : t70 = 0.1e1_dp/Qa
1297 634 : t73 = ba*x0a
1298 634 : t74 = t38 - x0a
1299 634 : t75 = t74**0.20e1_dp
1300 634 : t77 = LOG(t75*t61)
1301 634 : t80 = 0.20e1_dp*ba + 0.40e1_dp*x0a
1302 634 : t84 = x0a**0.20e1_dp
1303 634 : t86 = 0.1e1_dp/(t84 + t73 + ca)
1304 : t89 = t64 + 0.20e1_dp*ba*t68*t70 - t73*(t77 + t80*t68*t70) &
1305 634 : *t86
1306 634 : t90 = aa*t89
1307 634 : t91 = 0.1e1_dp/d2f0
1308 634 : t92 = t90*t91
1309 634 : t93 = t57*t92
1310 634 : t96 = t19 + 0.9531842930e0_dp*bf*t20 + cf
1311 634 : t97 = 0.1e1_dp/t96
1312 634 : t100 = LOG(0.9085602964e0_dp*t18*t97)
1313 634 : t101 = t28 + bf
1314 634 : t104 = ATAN(Qf/t101)
1315 634 : t106 = 0.1e1_dp/Qf
1316 634 : t109 = bf*x0f
1317 634 : t110 = t38 - x0f
1318 634 : t111 = t110**0.20e1_dp
1319 634 : t113 = LOG(t111*t97)
1320 634 : t116 = 0.20e1_dp*bf + 0.40e1_dp*x0f
1321 634 : t120 = x0f**0.20e1_dp
1322 634 : t122 = 0.1e1_dp/(t120 + t109 + cf)
1323 : t125 = t100 + 0.20e1_dp*bf*t104*t106 - t109*(t113 + t116*t104 &
1324 634 : *t106)*t122
1325 634 : t126 = af*t125
1326 634 : t127 = t13*t126
1327 :
1328 634 : e_0(ip) = e_0(ip) + ((t55 + t93 + t127)*t2)*sc
1329 :
1330 634 : IF (order >= 1) THEN
1331 634 : t129 = t5**0.333333333e0_dp
1332 634 : t130 = t2**2
1333 634 : t131 = 0.1e1_dp/t130
1334 634 : t132 = t1*t131
1335 634 : t133 = t3 - t132
1336 634 : t136 = t8**0.333333333e0_dp
1337 634 : t137 = -t133
1338 634 : t140 = 0.2564881399e1_dp*t129*t133 + 0.2564881399e1_dp*t136*t137
1339 634 : t141 = t140*t12
1340 634 : t142 = t4**0.30e1_dp
1341 634 : t143 = t11*t142
1342 634 : t144 = t143*t133
1343 634 : t147 = (-t141 - 0.40e1_dp*t144)*ap
1344 634 : t148 = t147*t54
1345 634 : t149 = t17**(-0.6666666668e0_dp)
1346 634 : t150 = t149*t24
1347 634 : t151 = t16*t131
1348 634 : t154 = t23**2
1349 634 : t155 = 0.1e1_dp/t154
1350 634 : t156 = t18*t155
1351 634 : t157 = t149*t16
1352 634 : t159 = 0.3028534320e0_dp*t157*t131
1353 634 : t160 = t17**(-0.8333333334e0_dp)
1354 634 : t161 = bp*t160
1355 634 : t164 = -t159 - 0.1588640488e0_dp*t161*t151
1356 634 : t167 = -0.3028534320e0_dp*t150*t151 - 0.9085602964e0_dp*t156*t164
1357 634 : t168 = t17**(-0.3333333332e0_dp)
1358 634 : t169 = t167*t168
1359 634 : t172 = t29**2
1360 634 : t173 = 0.1e1_dp/t172
1361 634 : t174 = bp*t173
1362 634 : t175 = t174*t160
1363 634 : t176 = Qp**2
1364 634 : t178 = 0.1e1_dp + t176*t173
1365 634 : t179 = 0.1e1_dp/t178
1366 634 : t180 = t151*t179
1367 634 : t183 = t39**0.10e1_dp
1368 634 : t184 = t183*t24
1369 634 : t185 = t160*t16
1370 634 : t186 = t185*t131
1371 634 : t189 = t40*t155
1372 634 : t191 = -0.3177280976e0_dp*t184*t186 - t189*t164
1373 634 : t192 = t39**(-0.20e1_dp)
1374 634 : t193 = t191*t192
1375 634 : t195 = t45*t173
1376 634 : t196 = t195*t160
1377 : t202 = 0.1100642416e1_dp*t169*t23 + 0.6354561950e0_dp*t175*t180 - &
1378 634 : t37*(t193*t23 + 0.3177280975e0_dp*t196*t180)*t51
1379 634 : t203 = t15*t202
1380 634 : t204 = t140*t56
1381 634 : t205 = t204*t92
1382 634 : t206 = t144*t92
1383 634 : t207 = 0.40e1_dp*t206
1384 634 : t208 = t149*t61
1385 634 : t211 = t60**2
1386 634 : t212 = 0.1e1_dp/t211
1387 634 : t213 = t18*t212
1388 634 : t214 = ba*t160
1389 634 : t217 = -t159 - 0.1588640488e0_dp*t214*t151
1390 634 : t220 = -0.3028534320e0_dp*t208*t151 - 0.9085602964e0_dp*t213*t217
1391 634 : t221 = t220*t168
1392 634 : t224 = t65**2
1393 634 : t225 = 0.1e1_dp/t224
1394 634 : t226 = ba*t225
1395 634 : t227 = t226*t160
1396 634 : t228 = Qa**2
1397 634 : t230 = 0.1e1_dp + t228*t225
1398 634 : t231 = 0.1e1_dp/t230
1399 634 : t232 = t151*t231
1400 634 : t235 = t74**0.10e1_dp
1401 634 : t236 = t235*t61
1402 634 : t239 = t75*t212
1403 634 : t241 = -0.3177280976e0_dp*t236*t186 - t239*t217
1404 634 : t242 = t74**(-0.20e1_dp)
1405 634 : t243 = t241*t242
1406 634 : t245 = t80*t225
1407 634 : t246 = t245*t160
1408 : t252 = 0.1100642416e1_dp*t221*t60 + 0.6354561950e0_dp*t227*t232 - &
1409 634 : t73*(t243*t60 + 0.3177280975e0_dp*t246*t232)*t86
1410 634 : t253 = aa*t252
1411 634 : t254 = t253*t91
1412 634 : t255 = t57*t254
1413 634 : t256 = t141*t126
1414 634 : t257 = t126*t133
1415 634 : t258 = t143*t257
1416 634 : t259 = 0.40e1_dp*t258
1417 634 : t260 = t149*t97
1418 634 : t263 = t96**2
1419 634 : t264 = 0.1e1_dp/t263
1420 634 : t265 = t18*t264
1421 634 : t266 = bf*t160
1422 634 : t269 = -t159 - 0.1588640488e0_dp*t266*t151
1423 634 : t272 = -0.3028534320e0_dp*t260*t151 - 0.9085602964e0_dp*t265*t269
1424 634 : t273 = t272*t168
1425 634 : t276 = t101**2
1426 634 : t277 = 0.1e1_dp/t276
1427 634 : t278 = bf*t277
1428 634 : t279 = t278*t160
1429 634 : t280 = Qf**2
1430 634 : t282 = 0.1e1_dp + t280*t277
1431 634 : t283 = 0.1e1_dp/t282
1432 634 : t284 = t151*t283
1433 634 : t287 = t110**0.10e1_dp
1434 634 : t288 = t287*t97
1435 634 : t291 = t111*t264
1436 634 : t293 = -0.3177280976e0_dp*t288*t186 - t291*t269
1437 634 : t294 = t110**(-0.20e1_dp)
1438 634 : t295 = t293*t294
1439 634 : t297 = t116*t277
1440 634 : t298 = t297*t160
1441 : t304 = 0.1100642416e1_dp*t273*t96 + 0.6354561950e0_dp*t279*t284 - &
1442 634 : t109*(t295*t96 + 0.3177280975e0_dp*t298*t284)*t122
1443 634 : t305 = af*t304
1444 634 : t306 = t13*t305
1445 :
1446 : e_a(ip) = e_a(ip) + ((t148 + t203 + t205 - t207 + t255 + t256 + t259 + t306)*t2 + &
1447 634 : t55 + t93 + t127)*sc
1448 :
1449 634 : t309 = -t3 - t132
1450 634 : t312 = -t309
1451 634 : t315 = 0.2564881399e1_dp*t129*t309 + 0.2564881399e1_dp*t136*t312
1452 634 : t316 = t315*t12
1453 634 : t317 = t143*t309
1454 634 : t320 = (-t316 - 0.40e1_dp*t317)*ap
1455 634 : t321 = t320*t54
1456 634 : t322 = t315*t56
1457 634 : t323 = t322*t92
1458 634 : t324 = t317*t92
1459 634 : t325 = 0.40e1_dp*t324
1460 634 : t326 = t316*t126
1461 634 : t327 = t126*t309
1462 634 : t328 = t143*t327
1463 634 : t329 = 0.40e1_dp*t328
1464 :
1465 : e_b(ip) = e_b(ip) + ((t321 + t203 + t323 - t325 + t255 + t326 + t329 + t306)*t2 + &
1466 634 : t55 + t93 + t127)*sc
1467 :
1468 : END IF
1469 634 : IF (order >= 2) THEN
1470 0 : t332 = t5**(-0.666666667e0_dp)
1471 0 : t333 = t133**2
1472 0 : t337 = 0.1e1_dp/t130/t2
1473 0 : t338 = t1*t337
1474 0 : t340 = -0.2e1_dp*t131 + 0.2e1_dp*t338
1475 0 : t343 = t8**(-0.666666667e0_dp)
1476 0 : t344 = t137**2
1477 0 : t347 = -t340
1478 : t350 = 0.8549604655e0_dp*t332*t333 + 0.2564881399e1_dp*t129*t340 &
1479 0 : + 0.8549604655e0_dp*t343*t344 + 0.2564881399e1_dp*t136*t347
1480 0 : t351 = t350*t12
1481 0 : t352 = t140*t142
1482 0 : t353 = t352*t133
1483 0 : t355 = t4**0.20e1_dp
1484 0 : t356 = t11*t355
1485 0 : t357 = t356*t333
1486 0 : t359 = t143*t340
1487 : t362 = (-t351 - 0.80e1_dp*t353 - 0.1200e2_dp*t357 - 0.40e1_dp*t359)* &
1488 0 : ap
1489 0 : t363 = t362*t54
1490 0 : t364 = t147*t202
1491 0 : t365 = 0.2e1_dp*t364
1492 0 : t366 = t17**(-0.1666666667e1_dp)
1493 0 : t367 = t366*t24
1494 0 : t368 = 0.3141592654e1_dp**2
1495 0 : t369 = 0.1e1_dp/t368
1496 0 : t370 = t130**2
1497 0 : t371 = 0.1e1_dp/t370
1498 0 : t372 = t369*t371
1499 0 : t373 = t367*t372
1500 0 : t375 = t149*t155
1501 0 : t376 = t151*t164
1502 0 : t379 = t16*t337
1503 0 : t383 = 0.1e1_dp/t154/t23
1504 0 : t384 = t18*t383
1505 0 : t385 = t164**2
1506 0 : t388 = t366*t369
1507 0 : t389 = t388*t371
1508 0 : t390 = 0.2019022880e0_dp*t389
1509 0 : t392 = 0.6057068640e0_dp*t157*t337
1510 0 : t393 = t17**(-0.1833333333e1_dp)
1511 0 : t394 = bp*t393
1512 : t399 = -t390 + t392 - 0.1323867073e0_dp*t394*t372 + 0.3177280976e0_dp &
1513 0 : *t161*t379
1514 : t402 = -0.2019022880e0_dp*t373 + 0.6057068640e0_dp*t375*t376 + 0.6057068640e0_dp &
1515 : *t150*t379 + 0.1817120593e1_dp*t384*t385 - 0.9085602964e0_dp &
1516 0 : *t156*t399
1517 0 : t403 = t402*t168
1518 0 : t406 = t17**(-0.1333333333e1_dp)
1519 0 : t407 = t167*t406
1520 0 : t408 = t23*t16
1521 0 : t409 = t408*t131
1522 0 : t415 = 0.1e1_dp/t172/t29
1523 0 : t416 = bp*t415
1524 0 : t417 = t416*t366
1525 0 : t418 = t372*t179
1526 0 : t421 = t174*t393
1527 0 : t424 = t379*t179
1528 0 : t427 = t172**2
1529 0 : t429 = 0.1e1_dp/t427/t29
1530 0 : t430 = bp*t429
1531 0 : t431 = t430*t366
1532 0 : t432 = t178**2
1533 0 : t434 = 0.1e1_dp/t432*t176
1534 0 : t435 = t372*t434
1535 0 : t439 = t183*t155
1536 0 : t440 = t439*t160
1537 0 : t443 = t393*t369
1538 0 : t444 = t443*t371
1539 0 : t447 = t185*t337
1540 0 : t450 = t40*t383
1541 : t454 = 0.5047557200e-1_dp*t373 + 0.6354561952e0_dp*t440*t376 - 0.2647734147e0_dp &
1542 : *t184*t444 + 0.6354561952e0_dp*t184*t447 + 0.2e1_dp* &
1543 0 : t450*t385 - t189*t399
1544 0 : t455 = t454*t192
1545 0 : t457 = t39**(-0.30e1_dp)
1546 0 : t458 = t191*t457
1547 0 : t459 = t458*t23
1548 0 : t463 = t45*t415
1549 0 : t464 = t463*t366
1550 0 : t467 = t195*t393
1551 0 : t472 = t45*t429
1552 0 : t473 = t472*t366
1553 : t479 = 0.1100642416e1_dp*t403*t23 + 0.3668808052e0_dp*t407*t409 + &
1554 : 0.1100642416e1_dp*t169*t164 + 0.4038045758e0_dp*t417*t418 + 0.5295468292e0_dp &
1555 : *t421*t418 - 0.1270912390e1_dp*t175*t424 - 0.4038045758e0_dp &
1556 : *t431*t435 - t37*(t455*t23 + 0.3177280976e0_dp*t459 &
1557 : *t186 + t193*t164 + 0.2019022879e0_dp*t464*t418 + 0.2647734146e0_dp &
1558 : *t467*t418 - 0.6354561950e0_dp*t196*t424 - 0.2019022879e0_dp* &
1559 0 : t473*t435)*t51
1560 0 : t480 = t15*t479
1561 0 : t481 = t350*t56
1562 0 : t482 = t481*t92
1563 0 : t483 = t353*t92
1564 0 : t484 = 0.80e1_dp*t483
1565 0 : t485 = t204*t254
1566 0 : t486 = 0.2e1_dp*t485
1567 0 : t487 = t357*t92
1568 0 : t488 = 0.1200e2_dp*t487
1569 0 : t489 = t359*t92
1570 0 : t490 = 0.40e1_dp*t489
1571 0 : t491 = t144*t254
1572 0 : t492 = 0.80e1_dp*t491
1573 0 : t493 = t366*t61
1574 0 : t494 = t493*t372
1575 0 : t496 = t149*t212
1576 0 : t497 = t151*t217
1577 0 : t503 = 0.1e1_dp/t211/t60
1578 0 : t504 = t18*t503
1579 0 : t505 = t217**2
1580 0 : t508 = ba*t393
1581 : t513 = -t390 + t392 - 0.1323867073e0_dp*t508*t372 + 0.3177280976e0_dp &
1582 0 : *t214*t379
1583 : t516 = -0.2019022880e0_dp*t494 + 0.6057068640e0_dp*t496*t497 + 0.6057068640e0_dp &
1584 : *t208*t379 + 0.1817120593e1_dp*t504*t505 - 0.9085602964e0_dp &
1585 0 : *t213*t513
1586 0 : t517 = t516*t168
1587 0 : t520 = t220*t406
1588 0 : t521 = t60*t16
1589 0 : t522 = t521*t131
1590 0 : t528 = 0.1e1_dp/t224/t65
1591 0 : t529 = ba*t528
1592 0 : t530 = t529*t366
1593 0 : t531 = t372*t231
1594 0 : t534 = t226*t393
1595 0 : t537 = t379*t231
1596 0 : t540 = t224**2
1597 0 : t542 = 0.1e1_dp/t540/t65
1598 0 : t543 = ba*t542
1599 0 : t544 = t543*t366
1600 0 : t545 = t230**2
1601 0 : t547 = 0.1e1_dp/t545*t228
1602 0 : t548 = t372*t547
1603 0 : t552 = t235*t212
1604 0 : t553 = t552*t160
1605 0 : t560 = t75*t503
1606 : t564 = 0.5047557200e-1_dp*t494 + 0.6354561952e0_dp*t553*t497 - 0.2647734147e0_dp &
1607 : *t236*t444 + 0.6354561952e0_dp*t236*t447 + 0.2e1_dp* &
1608 0 : t560*t505 - t239*t513
1609 0 : t565 = t564*t242
1610 0 : t567 = t74**(-0.30e1_dp)
1611 0 : t568 = t241*t567
1612 0 : t569 = t568*t60
1613 0 : t573 = t80*t528
1614 0 : t574 = t573*t366
1615 0 : t577 = t245*t393
1616 0 : t582 = t80*t542
1617 0 : t583 = t582*t366
1618 : t591 = aa*(0.1100642416e1_dp*t517*t60 + 0.3668808052e0_dp*t520* &
1619 : t522 + 0.1100642416e1_dp*t221*t217 + 0.4038045758e0_dp*t530*t531 &
1620 : + 0.5295468292e0_dp*t534*t531 - 0.1270912390e1_dp*t227*t537 - 0.4038045758e0_dp &
1621 : *t544*t548 - t73*(t565*t60 + 0.3177280976e0_dp* &
1622 : t569*t186 + t243*t217 + 0.2019022879e0_dp*t574*t531 + 0.2647734146e0_dp &
1623 : *t577*t531 - 0.6354561950e0_dp*t246*t537 - 0.2019022879e0_dp &
1624 0 : *t583*t548)*t86)*t91
1625 0 : t592 = t57*t591
1626 0 : t593 = t351*t126
1627 0 : t594 = t352*t257
1628 0 : t595 = 0.80e1_dp*t594
1629 0 : t596 = t141*t305
1630 0 : t597 = 0.2e1_dp*t596
1631 0 : t598 = t126*t333
1632 0 : t599 = t356*t598
1633 0 : t600 = 0.1200e2_dp*t599
1634 0 : t601 = t305*t133
1635 0 : t602 = t143*t601
1636 0 : t603 = 0.80e1_dp*t602
1637 0 : t604 = t126*t340
1638 0 : t605 = t143*t604
1639 0 : t606 = 0.40e1_dp*t605
1640 0 : t607 = t366*t97
1641 0 : t608 = t607*t372
1642 0 : t610 = t149*t264
1643 0 : t611 = t151*t269
1644 0 : t617 = 0.1e1_dp/t263/t96
1645 0 : t618 = t18*t617
1646 0 : t619 = t269**2
1647 0 : t622 = bf*t393
1648 : t627 = -t390 + t392 - 0.1323867073e0_dp*t622*t372 + 0.3177280976e0_dp &
1649 0 : *t266*t379
1650 : t630 = -0.2019022880e0_dp*t608 + 0.6057068640e0_dp*t610*t611 + 0.6057068640e0_dp &
1651 : *t260*t379 + 0.1817120593e1_dp*t618*t619 - 0.9085602964e0_dp &
1652 0 : *t265*t627
1653 0 : t631 = t630*t168
1654 0 : t634 = t272*t406
1655 0 : t635 = t96*t16
1656 0 : t636 = t635*t131
1657 0 : t642 = 0.1e1_dp/t276/t101
1658 0 : t643 = bf*t642
1659 0 : t644 = t643*t366
1660 0 : t645 = t372*t283
1661 0 : t648 = t278*t393
1662 0 : t651 = t379*t283
1663 0 : t654 = t276**2
1664 0 : t656 = 0.1e1_dp/t654/t101
1665 0 : t657 = bf*t656
1666 0 : t658 = t657*t366
1667 0 : t659 = t282**2
1668 0 : t661 = 0.1e1_dp/t659*t280
1669 0 : t662 = t372*t661
1670 0 : t666 = t287*t264
1671 0 : t667 = t666*t160
1672 0 : t674 = t111*t617
1673 : t678 = 0.5047557200e-1_dp*t608 + 0.6354561952e0_dp*t667*t611 - 0.2647734147e0_dp &
1674 : *t288*t444 + 0.6354561952e0_dp*t288*t447 + 0.2e1_dp* &
1675 0 : t674*t619 - t291*t627
1676 0 : t679 = t678*t294
1677 0 : t681 = t110**(-0.30e1_dp)
1678 0 : t682 = t293*t681
1679 0 : t683 = t682*t96
1680 0 : t687 = t116*t642
1681 0 : t688 = t687*t366
1682 0 : t691 = t297*t393
1683 0 : t696 = t116*t656
1684 0 : t697 = t696*t366
1685 : t704 = af*(0.1100642416e1_dp*t631*t96 + 0.3668808052e0_dp*t634* &
1686 : t636 + 0.1100642416e1_dp*t273*t269 + 0.4038045758e0_dp*t644*t645 &
1687 : + 0.5295468292e0_dp*t648*t645 - 0.1270912390e1_dp*t279*t651 - 0.4038045758e0_dp &
1688 : *t658*t662 - t109*(t679*t96 + 0.3177280976e0_dp &
1689 : *t683*t186 + t295*t269 + 0.2019022879e0_dp*t688*t645 + 0.2647734146e0_dp &
1690 : *t691*t645 - 0.6354561950e0_dp*t298*t651 - 0.2019022879e0_dp &
1691 0 : *t697*t662)*t122)
1692 0 : t705 = t13*t704
1693 : t706 = t363 + t365 + t480 + t482 - t484 + t486 - t488 - t490 - t492 &
1694 0 : + t592 + t593 + t595 + t597 + t600 + t603 + t606 + t705
1695 0 : t709 = 0.2e1_dp*t203
1696 0 : t712 = 0.2e1_dp*t255
1697 0 : t715 = 0.2e1_dp*t306
1698 :
1699 : e_aa(ip) = e_aa(ip) + (t706*t2 + 0.2e1_dp*t148 + t709 + 0.2e1_dp*t205 - 0.80e1_dp*t206 &
1700 0 : + t712 + 0.2e1_dp*t256 + 0.80e1_dp*t258 + t715)*sc
1701 :
1702 0 : t716 = t332*t133
1703 0 : t719 = t129*t1
1704 0 : t722 = t343*t137
1705 0 : t725 = t136*t1
1706 : t728 = 0.8549604655e0_dp*t716*t309 + 0.5129762798e1_dp*t719*t337 &
1707 0 : + 0.8549604655e0_dp*t722*t312 - 0.5129762798e1_dp*t725*t337
1708 0 : t729 = t728*t12
1709 0 : t730 = t352*t309
1710 0 : t732 = t315*t142
1711 0 : t733 = t732*t133
1712 0 : t735 = t133*t309
1713 : t741 = (-t729 - 0.40e1_dp*t730 - 0.40e1_dp*t733 - 0.1200e2_dp*t356*t735 &
1714 0 : - 0.80e1_dp*t143*t338)*ap
1715 0 : t742 = t741*t54
1716 0 : t743 = t320*t202
1717 0 : t744 = t728*t56
1718 0 : t745 = t744*t92
1719 0 : t746 = t730*t92
1720 0 : t748 = t733*t92
1721 0 : t750 = t356*t133
1722 0 : t751 = t91*t309
1723 0 : t752 = t90*t751
1724 0 : t753 = t750*t752
1725 0 : t755 = t143*t1
1726 0 : t756 = t337*aa
1727 0 : t757 = t89*t91
1728 0 : t758 = t756*t757
1729 0 : t759 = t755*t758
1730 0 : t762 = t322*t254
1731 : t763 = t742 + t364 + t743 + t480 + t745 - 0.40e1_dp*t746 + t485 - 0.40e1_dp &
1732 0 : *t748 - 0.1200e2_dp*t753 - 0.80e1_dp*t759 - 0.40e1_dp*t491 + t762
1733 0 : t764 = t317*t254
1734 0 : t766 = t729*t126
1735 0 : t767 = t352*t327
1736 0 : t769 = t732*t257
1737 0 : t771 = t356*af
1738 0 : t772 = t125*t133
1739 0 : t773 = t772*t309
1740 0 : t774 = t771*t773
1741 0 : t777 = t143*af
1742 0 : t778 = t125*t1
1743 0 : t779 = t778*t337
1744 0 : t780 = t777*t779
1745 0 : t782 = t316*t305
1746 0 : t783 = t305*t309
1747 0 : t784 = t143*t783
1748 : t786 = -0.40e1_dp*t764 + t592 + t766 + 0.40e1_dp*t767 + t596 + 0.40e1_dp &
1749 : *t769 + 0.1200e2_dp*t774 + 0.40e1_dp*t602 + 0.80e1_dp*t780 + t782 + &
1750 0 : 0.40e1_dp*t784 + t705
1751 :
1752 : e_ab(ip) = e_ab(ip) + ((t763 + t786)*t2 + t148 + t709 + t205 - t207 + t712 + t256 &
1753 0 : + t259 + t715 + t321 + t323 - t325 + t326 + t329)*sc
1754 0 : t789 = t309**2
1755 0 : t793 = 0.2e1_dp*t131 + 0.2e1_dp*t338
1756 0 : t796 = t312**2
1757 0 : t799 = -t793
1758 : t802 = 0.8549604655e0_dp*t332*t789 + 0.2564881399e1_dp*t129*t793 &
1759 0 : + 0.8549604655e0_dp*t343*t796 + 0.2564881399e1_dp*t136*t799
1760 0 : t803 = t802*t12
1761 0 : t804 = t732*t309
1762 0 : t806 = t356*t789
1763 0 : t808 = t143*t793
1764 : t811 = (-t803 - 0.80e1_dp*t804 - 0.1200e2_dp*t806 - 0.40e1_dp*t808)* &
1765 0 : ap
1766 0 : t812 = t811*t54
1767 0 : t813 = 0.2e1_dp*t743
1768 0 : t814 = t802*t56
1769 0 : t815 = t814*t92
1770 0 : t816 = t804*t92
1771 0 : t817 = 0.80e1_dp*t816
1772 0 : t818 = 0.2e1_dp*t762
1773 0 : t819 = t806*t92
1774 0 : t820 = 0.1200e2_dp*t819
1775 0 : t821 = t808*t92
1776 0 : t822 = 0.40e1_dp*t821
1777 0 : t823 = 0.80e1_dp*t764
1778 0 : t824 = t803*t126
1779 0 : t825 = t732*t327
1780 0 : t826 = 0.80e1_dp*t825
1781 0 : t827 = 0.2e1_dp*t782
1782 0 : t828 = t126*t789
1783 0 : t829 = t356*t828
1784 0 : t830 = 0.1200e2_dp*t829
1785 0 : t831 = 0.80e1_dp*t784
1786 0 : t832 = t126*t793
1787 0 : t833 = t143*t832
1788 0 : t834 = 0.40e1_dp*t833
1789 : t835 = t812 + t813 + t480 + t815 - t817 + t818 - t820 - t822 - t823 &
1790 0 : + t592 + t824 + t826 + t827 + t830 + t831 + t834 + t705
1791 :
1792 : e_bb(ip) = e_bb(ip) + (t835*t2 + 0.2e1_dp*t321 + t709 + 0.2e1_dp*t323 - 0.80e1_dp*t324 &
1793 0 : + t712 + 0.2e1_dp*t326 + 0.80e1_dp*t328 + t715)*sc
1794 :
1795 : END IF
1796 634 : IF (order >= 3) THEN
1797 0 : t842 = 0.3e1_dp*t480
1798 0 : t851 = 0.3e1_dp*t592
1799 0 : t857 = 0.3e1_dp*t705
1800 0 : t859 = t17**(-0.2666666667e1_dp)
1801 0 : t862 = 0.1e1_dp/t368/0.3141592654e1_dp
1802 0 : t864 = 0.1e1_dp/t370/t130
1803 0 : t865 = t862*t864
1804 0 : t866 = t859*t24*t865
1805 0 : t869 = t372*t164
1806 0 : t870 = t366*t155*t869
1807 0 : t873 = 0.1e1_dp/t370/t2
1808 0 : t874 = t369*t873
1809 0 : t875 = t367*t874
1810 0 : t878 = t151*t385
1811 0 : t881 = t379*t164
1812 0 : t884 = t151*t399
1813 0 : t887 = t16*t371
1814 0 : t890 = t154**2
1815 0 : t891 = 0.1e1_dp/t890
1816 0 : t893 = t385*t164
1817 0 : t896 = t164*t399
1818 0 : t901 = 0.3365038134e0_dp*t859*t862*t864
1819 0 : t903 = 0.1211413728e1_dp*t388*t873
1820 0 : t905 = 0.1817120592e1_dp*t157*t371
1821 0 : t906 = t17**(-0.2833333333e1_dp)
1822 : t914 = -t901 + t903 - t905 - 0.2427089633e0_dp*bp*t906*t865 + 0.7943202439e0_dp &
1823 0 : *t394*t874 - 0.9531842928e0_dp*t161*t887
1824 0 : t937 = t17**(-0.2666666666e1_dp)
1825 0 : t942 = t906*t862*t864
1826 0 : t945 = t443*t873
1827 0 : t948 = t185*t371
1828 : t957 = 0.8412595335e-1_dp*t866 - 0.1514267160e0_dp*t870 - 0.3028534320e0_dp &
1829 : *t875 - 0.1906368585e1_dp*t183*t383*t160*t878 + 0.7943202441e0_dp &
1830 : *t439*t393*t869 - 0.1906368585e1_dp*t440*t881 + 0.9531842928e0_dp &
1831 : *t440*t884 + 0.4206297667e-1_dp*t24*t937*t865 - 0.4854179269e0_dp &
1832 : *t184*t942 + 0.1588640488e1_dp*t184*t945 - 0.1906368586e1_dp &
1833 : *t184*t948 - 0.6e1_dp*t40*t891*t893 + 0.6e1_dp*t450 &
1834 0 : *t896 - t189*t914
1835 0 : t972 = t39**(-0.40e1_dp)
1836 0 : t979 = t874*t179
1837 0 : t982 = 0.1e1_dp/t427
1838 0 : t984 = t17**(-0.2500000000e1_dp)
1839 0 : t986 = t865*t179
1840 0 : t993 = 0.1e1_dp/t427/t172
1841 0 : t996 = t865*t434
1842 0 : t1010 = t887*t179
1843 0 : t1013 = t874*t434
1844 0 : t1019 = t427**2
1845 0 : t1020 = 0.1e1_dp/t1019
1846 0 : t1025 = t176**2
1847 0 : t1027 = t865/t432/t178*t1025
1848 : t1030 = t957*t192*t23 + 0.2e1_dp*t455*t164 + 0.6354561952e0_dp* &
1849 : t454*t457*t23*t186 + t193*t399 + 0.6354561952e0_dp*t458*t164 &
1850 : *t186 - 0.6354561952e0_dp*t459*t447 + 0.1514267160e0_dp*t191 &
1851 : *t972*t23*t389 + 0.2647734147e0_dp*t459*t444 - 0.1211413727e1_dp &
1852 : *t464*t979 + 0.1924500894e0_dp*t45*t982*t984*t986 + 0.3365038132e0_dp &
1853 : *t463*t859*t986 - 0.4490502088e0_dp*t45*t993*t984 &
1854 : *t996 - 0.1588640487e1_dp*t467*t979 + 0.1682519066e0_dp*t463* &
1855 : t937*t986 + 0.4854179267e0_dp*t195*t906*t986 - 0.1682519066e0_dp &
1856 : *t472*t937*t996 + 0.1906368585e1_dp*t196*t1010 + 0.1211413727e1_dp &
1857 : *t473*t1013 - 0.3365038132e0_dp*t472*t859*t996 + 0.2566001193e0_dp &
1858 0 : *t45*t1020*t984*t1027
1859 0 : t1043 = t17**(-0.2333333333e1_dp)
1860 : t1084 = 0.1100642416e1_dp*(-0.3365038134e0_dp*t866 + 0.6057068641e0_dp* &
1861 : t870 + 0.1211413728e1_dp*t875 - 0.1817120592e1_dp*t149*t383*t878 &
1862 : - 0.1817120592e1_dp*t375*t881 + 0.9085602960e0_dp*t375*t884 - &
1863 : 0.1817120592e1_dp*t150*t887 - 0.5451361779e1_dp*t18*t891*t893 &
1864 : + 0.5451361779e1_dp*t384*t896 - 0.9085602964e0_dp*t156*t914)*t168 &
1865 : *t23 + 0.2201284832e1_dp*t403*t164 - t37*t1030*t51 + 0.7337616104e0_dp &
1866 : *t402*t406*t409 + 0.1100642416e1_dp*t169*t399 + &
1867 : 0.7337616104e0_dp*t407*t376 - 0.7337616104e0_dp*t407*t408*t337 &
1868 : + 0.4891744068e0_dp*t167*t1043*t23*t369*t371 - 0.2422827454e1_dp &
1869 : *t417*t979 + 0.3849001789e0_dp*bp*t982*t984*t986 + 0.6730076265e0_dp &
1870 : *t416*t859*t986 - 0.8981004177e0_dp*bp*t993*t984 &
1871 : *t996 - 0.3177280975e1_dp*t421*t979 + 0.3365038132e0_dp*t416* &
1872 : t937*t986 + 0.9708358534e0_dp*t174*t906*t986 - 0.3365038132e0_dp &
1873 : *t430*t937*t996 + 0.3812737170e1_dp*t175*t1010 + 0.2422827454e1_dp &
1874 : *t431*t1013 - 0.6730076265e0_dp*t430*t859*t996 + 0.5132002385e0_dp &
1875 0 : *bp*t1020*t984*t1027
1876 0 : t1085 = t15*t1084
1877 0 : t1086 = t147*t479
1878 0 : t1088 = t5**(-0.1666666667e1_dp)
1879 0 : t1089 = t333*t133
1880 0 : t1094 = t1*t371
1881 0 : t1096 = 0.6e1_dp*t337 - 0.6e1_dp*t1094
1882 0 : t1099 = t8**(-0.1666666667e1_dp)
1883 : t1108 = -0.5699736440e0_dp*t1088*t1089 + 0.2564881396e1_dp*t716*t340 &
1884 : + 0.2564881399e1_dp*t129*t1096 - 0.5699736440e0_dp*t1099*t344 &
1885 : *t137 + 0.2564881396e1_dp*t722*t347 - 0.2564881399e1_dp*t136* &
1886 0 : t1096
1887 0 : t1109 = t1108*t12
1888 0 : t1110 = t350*t142
1889 0 : t1111 = t1110*t133
1890 0 : t1113 = t140*t355
1891 0 : t1114 = t1113*t333
1892 0 : t1116 = t352*t340
1893 0 : t1118 = t4**0.10e1_dp
1894 0 : t1119 = t11*t1118
1895 0 : t1120 = t1119*t1089
1896 0 : t1125 = t143*t1096
1897 0 : t1137 = t351*t305
1898 0 : t1142 = t352*t601
1899 0 : t1145 = t859*t97*t865
1900 0 : t1148 = t372*t269
1901 0 : t1149 = t366*t264*t1148
1902 0 : t1151 = t607*t874
1903 0 : t1154 = t151*t619
1904 0 : t1157 = t379*t269
1905 0 : t1160 = t151*t627
1906 0 : t1165 = t263**2
1907 0 : t1166 = 0.1e1_dp/t1165
1908 0 : t1168 = t619*t269
1909 0 : t1171 = t269*t627
1910 : t1181 = -t901 + t903 - t905 - 0.2427089633e0_dp*bf*t906*t865 + 0.7943202439e0_dp &
1911 0 : *t622*t874 - 0.9531842928e0_dp*t266*t887
1912 0 : t1205 = t874*t283
1913 0 : t1208 = 0.1e1_dp/t654
1914 0 : t1211 = t865*t283
1915 0 : t1218 = 0.1e1_dp/t654/t276
1916 0 : t1221 = t865*t661
1917 0 : t1235 = t887*t283
1918 0 : t1238 = t874*t661
1919 0 : t1244 = t654**2
1920 0 : t1245 = 0.1e1_dp/t1244
1921 0 : t1250 = t280**2
1922 0 : t1252 = t865/t659/t282*t1250
1923 : t1284 = 0.8412595335e-1_dp*t1145 - 0.1514267160e0_dp*t1149 - 0.3028534320e0_dp &
1924 : *t1151 - 0.1906368585e1_dp*t287*t617*t160*t1154 + 0.7943202441e0_dp &
1925 : *t666*t393*t1148 - 0.1906368585e1_dp*t667*t1157 &
1926 : + 0.9531842928e0_dp*t667*t1160 + 0.4206297667e-1_dp*t97*t937*t865 &
1927 : - 0.4854179269e0_dp*t288*t942 + 0.1588640488e1_dp*t288*t945 &
1928 : - 0.1906368586e1_dp*t288*t948 - 0.6e1_dp*t111*t1166*t1168 + 0.6e1_dp &
1929 0 : *t674*t1171 - t291*t1181
1930 0 : t1299 = t110**(-0.40e1_dp)
1931 : t1341 = t1284*t294*t96 + 0.2e1_dp*t679*t269 + 0.6354561952e0_dp* &
1932 : t678*t681*t96*t186 + t295*t627 + 0.6354561952e0_dp*t682* &
1933 : t269*t186 - 0.6354561952e0_dp*t683*t447 + 0.1514267160e0_dp*t293 &
1934 : *t1299*t96*t389 + 0.2647734147e0_dp*t683*t444 - 0.1211413727e1_dp &
1935 : *t688*t1205 + 0.1924500894e0_dp*t116*t1208*t984*t1211 &
1936 : + 0.3365038132e0_dp*t687*t859*t1211 - 0.4490502088e0_dp*t116*t1218 &
1937 : *t984*t1221 - 0.1588640487e1_dp*t691*t1205 + 0.1682519066e0_dp &
1938 : *t687*t937*t1211 + 0.4854179267e0_dp*t297*t906*t1211 - &
1939 : 0.1682519066e0_dp*t696*t937*t1221 + 0.1906368585e1_dp*t298*t1235 &
1940 : + 0.1211413727e1_dp*t697*t1238 - 0.3365038132e0_dp*t696*t859 &
1941 0 : *t1221 + 0.2566001193e0_dp*t116*t1245*t984*t1252
1942 : t1344 = 0.1100642416e1_dp*(-0.3365038134e0_dp*t1145 + 0.6057068641e0_dp &
1943 : *t1149 + 0.1211413728e1_dp*t1151 - 0.1817120592e1_dp*t149*t617* &
1944 : t1154 - 0.1817120592e1_dp*t610*t1157 + 0.9085602960e0_dp*t610*t1160 &
1945 : - 0.1817120592e1_dp*t260*t887 - 0.5451361779e1_dp*t18*t1166 &
1946 : *t1168 + 0.5451361779e1_dp*t618*t1171 - 0.9085602964e0_dp*t265* &
1947 : t1181)*t168*t96 + 0.2201284832e1_dp*t631*t269 + 0.7337616104e0_dp &
1948 : *t630*t406*t636 + 0.1100642416e1_dp*t273*t627 + 0.7337616104e0_dp &
1949 : *t634*t611 - 0.7337616104e0_dp*t634*t635*t337 + 0.4891744068e0_dp &
1950 : *t272*t1043*t96*t369*t371 - 0.2422827454e1_dp*t644 &
1951 : *t1205 + 0.3849001789e0_dp*bf*t1208*t984*t1211 + 0.6730076265e0_dp &
1952 : *t643*t859*t1211 - 0.8981004177e0_dp*bf*t1218*t984* &
1953 : t1221 - 0.3177280975e1_dp*t648*t1205 + 0.3365038132e0_dp*t643*t937 &
1954 : *t1211 + 0.9708358534e0_dp*t278*t906*t1211 - 0.3365038132e0_dp &
1955 : *t657*t937*t1221 + 0.3812737170e1_dp*t279*t1235 + 0.2422827454e1_dp &
1956 : *t658*t1238 - 0.6730076265e0_dp*t657*t859*t1221 + 0.5132002385e0_dp &
1957 0 : *bf*t1245*t984*t1252 - t109*t1341*t122
1958 0 : t1346 = t13*af*t1344
1959 0 : t1358 = t356*t305*t333
1960 : t1360 = t1085 + 0.3e1_dp*t1086 + (-t1109 - 0.120e2_dp*t1111 - 0.3600e2_dp &
1961 : *t1114 - 0.120e2_dp*t1116 - 0.24000e2_dp*t1120 - 0.3600e2_dp*t356 &
1962 : *t133*t340 - 0.40e1_dp*t1125)*ap*t54 + t1109*t126 - 0.24000e2_dp &
1963 : *t1120*t92 + 0.120e2_dp*t1110*t257 - 0.120e2_dp*t1116*t92 &
1964 : + 0.3e1_dp*t1137 + 0.3600e2_dp*t771*t772*t340 + 0.240e2_dp*t1142 &
1965 : + t1346 - 0.120e2_dp*t1111*t92 - 0.3600e2_dp*t1114*t92 - 0.3600e2_dp &
1966 0 : *t750*t90*t91*t340 - 0.40e1_dp*t1125*t92 + 0.3600e2_dp*t1358
1967 0 : t1361 = t362*t202
1968 0 : t1363 = t353*t254
1969 0 : t1365 = t144*t591
1970 0 : t1368 = t859*t61*t865
1971 0 : t1371 = t372*t217
1972 0 : t1372 = t366*t212*t1371
1973 0 : t1374 = t493*t874
1974 0 : t1377 = t151*t505
1975 0 : t1380 = t379*t217
1976 0 : t1383 = t151*t513
1977 0 : t1388 = t211**2
1978 0 : t1389 = 0.1e1_dp/t1388
1979 0 : t1391 = t505*t217
1980 0 : t1394 = t217*t513
1981 : t1404 = -t901 + t903 - t905 - 0.2427089633e0_dp*ba*t906*t865 + 0.7943202439e0_dp &
1982 0 : *t508*t874 - 0.9531842928e0_dp*t214*t887
1983 : t1455 = 0.8412595335e-1_dp*t1368 - 0.1514267160e0_dp*t1372 - 0.3028534320e0_dp &
1984 : *t1374 - 0.1906368585e1_dp*t235*t503*t160*t1377 + 0.7943202441e0_dp &
1985 : *t552*t393*t1371 - 0.1906368585e1_dp*t553*t1380 &
1986 : + 0.9531842928e0_dp*t553*t1383 + 0.4206297667e-1_dp*t61*t937*t865 &
1987 : - 0.4854179269e0_dp*t236*t942 + 0.1588640488e1_dp*t236*t945 &
1988 : - 0.1906368586e1_dp*t236*t948 - 0.6e1_dp*t75*t1389*t1391 + 0.6e1_dp &
1989 0 : *t560*t1394 - t239*t1404
1990 0 : t1469 = t74**(-0.40e1_dp)
1991 0 : t1477 = t874*t231
1992 0 : t1480 = 0.1e1_dp/t540
1993 0 : t1483 = t865*t231
1994 0 : t1490 = 0.1e1_dp/t540/t224
1995 0 : t1493 = t865*t547
1996 0 : t1507 = t887*t231
1997 0 : t1510 = t874*t547
1998 0 : t1516 = t540**2
1999 0 : t1517 = 0.1e1_dp/t1516
2000 0 : t1522 = t228**2
2001 0 : t1524 = t865/t545/t230*t1522
2002 : t1527 = t1455*t242*t60 + 0.2e1_dp*t565*t217 + 0.6354561952e0_dp* &
2003 : t564*t567*t60*t186 + 0.6354561952e0_dp*t568*t217*t186 - &
2004 : 0.6354561952e0_dp*t569*t447 + 0.1514267160e0_dp*t241*t1469*t60 &
2005 : *t389 + 0.2647734147e0_dp*t569*t444 + t243*t513 - 0.1211413727e1_dp &
2006 : *t574*t1477 + 0.1924500894e0_dp*t80*t1480*t984*t1483 + &
2007 : 0.3365038132e0_dp*t573*t859*t1483 - 0.4490502088e0_dp*t80*t1490 &
2008 : *t984*t1493 - 0.1588640487e1_dp*t577*t1477 + 0.1682519066e0_dp &
2009 : *t573*t937*t1483 + 0.4854179267e0_dp*t245*t906*t1483 - 0.1682519066e0_dp &
2010 : *t582*t937*t1493 + 0.1906368585e1_dp*t246*t1507 &
2011 : + 0.1211413727e1_dp*t583*t1510 - 0.3365038132e0_dp*t582*t859* &
2012 0 : t1493 + 0.2566001193e0_dp*t80*t1517*t984*t1524
2013 : t1567 = 0.1100642416e1_dp*(-0.3365038134e0_dp*t1368 + 0.6057068641e0_dp &
2014 : *t1372 + 0.1211413728e1_dp*t1374 - 0.1817120592e1_dp*t149*t503* &
2015 : t1377 - 0.1817120592e1_dp*t496*t1380 + 0.9085602960e0_dp*t496*t1383 &
2016 : - 0.1817120592e1_dp*t208*t887 - 0.5451361779e1_dp*t18*t1389 &
2017 : *t1391 + 0.5451361779e1_dp*t504*t1394 - 0.9085602964e0_dp*t213* &
2018 : t1404)*t168*t60 + 0.2201284832e1_dp*t517*t217 + 0.7337616104e0_dp &
2019 : *t516*t406*t522 + 0.7337616104e0_dp*t520*t497 - 0.7337616104e0_dp &
2020 : *t520*t521*t337 + 0.4891744068e0_dp*t220*t1043*t60* &
2021 : t369*t371 - t73*t1527*t86 + 0.1100642416e1_dp*t221*t513 - 0.2422827454e1_dp &
2022 : *t530*t1477 + 0.3849001789e0_dp*ba*t1480*t984 &
2023 : *t1483 + 0.6730076265e0_dp*t529*t859*t1483 - 0.8981004177e0_dp* &
2024 : ba*t1490*t984*t1493 - 0.3177280975e1_dp*t534*t1477 + 0.3365038132e0_dp &
2025 : *t529*t937*t1483 + 0.9708358534e0_dp*t226*t906*t1483 &
2026 : - 0.3365038132e0_dp*t543*t937*t1493 + 0.3812737170e1_dp*t227 &
2027 : *t1507 + 0.2422827454e1_dp*t544*t1510 - 0.6730076265e0_dp*t543* &
2028 0 : t859*t1493 + 0.5132002385e0_dp*ba*t1517*t984*t1524
2029 0 : t1570 = t57*aa*t1567*t91
2030 0 : t1576 = t481*t254
2031 0 : t1578 = t357*t254
2032 0 : t1580 = t204*t591
2033 0 : t1582 = t359*t254
2034 0 : t1588 = t143*t305*t340
2035 0 : t1590 = t141*t704
2036 0 : t1593 = t143*t704*t133
2037 : t1599 = 0.3e1_dp*t1361 - 0.240e2_dp*t1363 - 0.120e2_dp*t1365 + t1570 + &
2038 : 0.40e1_dp*t143*t126*t1096 + t1108*t56*t92 + 0.3e1_dp*t1576 &
2039 : - 0.3600e2_dp*t1578 + 0.3e1_dp*t1580 - 0.120e2_dp*t1582 + 0.24000e2_dp* &
2040 : t1119*t126*t1089 + 0.120e2_dp*t1588 + 0.3e1_dp*t1590 + 0.120e2_dp &
2041 0 : *t1593 + 0.120e2_dp*t352*t604 + 0.3600e2_dp*t1113*t598
2042 :
2043 : e_aaa(ip) = e_aaa(ip) + (t842 + 0.6e1_dp*t364 + 0.3e1_dp*t363 + 0.3e1_dp*t482 + 0.6e1_dp* &
2044 : t485 - 0.240e2_dp*t483 - 0.3600e2_dp*t487 - 0.120e2_dp*t489 - 0.240e2_dp &
2045 : *t491 + t851 + 0.3e1_dp*t593 + 0.6e1_dp*t596 + 0.240e2_dp*t594 + 0.3600e2_dp &
2046 : *t599 + 0.240e2_dp*t602 + t857 + 0.120e2_dp*t605 + (t1360 + &
2047 0 : t1599)*t2)*sc
2048 :
2049 0 : t1602 = 0.2400e2_dp*t753
2050 0 : t1603 = 0.2e1_dp*t766
2051 0 : t1604 = 0.80e1_dp*t746
2052 0 : t1605 = 0.2e1_dp*t745
2053 0 : t1606 = 0.80e1_dp*t748
2054 0 : t1609 = 0.160e2_dp*t759
2055 : t1610 = -t1602 + t1603 - t1604 + t1605 + t818 - t490 - t1606 + t827 &
2056 : + t363 - t488 + t831 + t813 - t823 - 0.160e2_dp*t491 + 0.4e1_dp*t364 &
2057 0 : + t842 - t1609
2058 0 : t1611 = 0.80e1_dp*t767
2059 0 : t1613 = t322*t591
2060 0 : t1615 = 0.80e1_dp*t732*t601
2061 0 : t1620 = 0.80e1_dp*t733*t254
2062 0 : t1622 = 0.80e1_dp*t352*t783
2063 0 : t1626 = t732*t340
2064 0 : t1639 = 0.2e1_dp*t337 - 0.6e1_dp*t1094
2065 : t1653 = -0.5699736440e0_dp*t1088*t333*t309 + 0.3419841862e1_dp*t716 &
2066 : *t338 + 0.8549604655e0_dp*t332*t340*t309 + 0.2564881399e1_dp* &
2067 : t129*t1639 - 0.5699736440e0_dp*t1099*t344*t312 - 0.3419841862e1_dp &
2068 : *t722*t338 + 0.8549604655e0_dp*t343*t347*t312 - 0.2564881399e1_dp &
2069 0 : *t136*t1639
2070 0 : t1654 = t1653*t12
2071 0 : t1657 = t143*t704*t309
2072 : t1659 = t1085 + 0.2e1_dp*t1086 + t1613 + t1615 - 0.160e2_dp*t352*t1 &
2073 : *t758 - t1620 + t1622 + 0.40e1_dp*t1110*t327 + t1137 + 0.80e1_dp* &
2074 0 : t1142 - 0.40e1_dp*t1626*t92 + t1346 + t1654*t126 + 0.40e1_dp*t1657
2075 0 : t1664 = 0.160e2_dp*t777*t304*t1*t337
2076 0 : t1666 = 0.80e1_dp*t730*t254
2077 0 : t1668 = t143*t1639
2078 0 : t1671 = t728*t142
2079 0 : t1672 = t1671*t133
2080 0 : t1675 = t1110*t309
2081 0 : t1690 = 0.2400e2_dp*t771*t304*t133*t309
2082 : t1694 = 0.1200e2_dp*t1358 + t1664 + t1361 - t1666 - 0.80e1_dp*t1363 - &
2083 : 0.40e1_dp*t1668*t92 - 0.80e1_dp*t1672*t92 - 0.40e1_dp*t1675*t92 &
2084 : - 0.2400e2_dp*t1113*t133*t752 - 0.80e1_dp*t1365 + t1570 - 0.4800e2_dp &
2085 : *t356*t133*aa*t757*t338 + t1690 + 0.40e1_dp*t143*t126 &
2086 0 : *t1639
2087 0 : t1696 = t316*t704
2088 0 : t1697 = t315*t355
2089 0 : t1698 = t1697*t333
2090 0 : t1701 = t1119*af
2091 : t1726 = t1696 - 0.1200e2_dp*t1698*t92 + 0.24000e2_dp*t1701*t125* &
2092 : t333*t309 + t1653*t56*t92 + 0.2400e2_dp*t1113*af*t773 + &
2093 : 0.4800e2_dp*t771*t772*t338 + 0.160e2_dp*t352*af*t779 + 0.40e1_dp &
2094 : *t732*t604 + t1576 - 0.1200e2_dp*t1578 + 0.1200e2_dp*t1697*t598 &
2095 0 : + 0.2e1_dp*t1580 - 0.40e1_dp*t1582 + 0.80e1_dp*t1671*t257
2096 0 : t1730 = 0.160e2_dp*t755*t756*t252*t91
2097 : t1750 = -t1654 - 0.40e1_dp*t1675 - 0.80e1_dp*t1672 - 0.2400e2_dp*t1113 &
2098 : *t735 - 0.160e2_dp*t352*t338 - 0.1200e2_dp*t1698 - 0.24000e2_dp*t1119 &
2099 : *t333*t309 - 0.4800e2_dp*t356*t133*t1*t337 - 0.40e1_dp* &
2100 0 : t1626 - 0.1200e2_dp*t356*t340*t309 - 0.40e1_dp*t1668
2101 0 : t1754 = 0.2e1_dp*t744*t254
2102 0 : t1764 = 0.2e1_dp*t741*t202
2103 0 : t1765 = t320*t479
2104 0 : t1768 = 0.2400e2_dp*t750*t253*t751
2105 0 : t1775 = 0.2e1_dp*t729*t305
2106 0 : t1776 = t317*t591
2107 : t1778 = -t1730 + t1750*ap*t54 + t1754 - 0.24000e2_dp*t1119*t333 &
2108 : *t752 + 0.40e1_dp*t1588 - 0.1200e2_dp*t356*t340*t752 + 0.2e1_dp &
2109 : *t1590 + t1764 + t1765 - t1768 + 0.80e1_dp*t1593 + 0.1200e2_dp*t771 &
2110 0 : *t125*t340*t309 + t1775 - 0.40e1_dp*t1776
2111 0 : t1782 = 0.2e1_dp*t742
2112 0 : t1783 = 0.160e2_dp*t780
2113 0 : t1785 = 0.2400e2_dp*t774
2114 0 : t1786 = 0.80e1_dp*t769
2115 : t1789 = t606 + t1611 + t600 + t482 + t851 + t595 + (t1659 + t1694 + &
2116 : t1726 + t1778)*t2 + t1782 + t1783 + t593 + 0.4e1_dp*t596 + t857 &
2117 0 : - t484 + t1785 + t1786 + 0.4e1_dp*t485 + 0.160e2_dp*t602
2118 :
2119 0 : e_aab(ip) = e_aab(ip) + (t1610 + t1789)*sc
2120 :
2121 : t1795 = -t1602 + t826 + t812 + t824 + t834 + t1603 - t1604 + t1605 &
2122 : + 0.4e1_dp*t762 - t1606 + 0.4e1_dp*t782 + t830 + 0.160e2_dp*t784 + 0.4e1_dp &
2123 0 : *t743 - 0.160e2_dp*t764 - t492 + t365
2124 0 : t1797 = t143*t305*t793
2125 0 : t1804 = t337*t309
2126 : t1826 = -0.5699736440e0_dp*t1088*t133*t789 + 0.3419841862e1_dp*t332 &
2127 : *t1*t1804 + 0.8549604655e0_dp*t716*t793 - 0.5129762798e1_dp* &
2128 : t129*t337 - 0.1538928839e2_dp*t719*t371 - 0.5699736440e0_dp*t1099 &
2129 : *t137*t796 - 0.3419841862e1_dp*t343*t1*t337*t312 + 0.8549604655e0_dp &
2130 : *t722*t799 + 0.5129762798e1_dp*t136*t337 + 0.1538928839e2_dp &
2131 0 : *t725*t371
2132 0 : t1829 = t352*t793
2133 0 : t1832 = t814*t254
2134 0 : t1836 = t732*t783
2135 0 : t1838 = t811*t202
2136 : t1839 = t1085 + t1086 + 0.40e1_dp*t1797 + 0.2e1_dp*t1613 + t1615 + t1826 &
2137 : *t56*t92 - 0.40e1_dp*t1829*t92 - t1620 + t1832 + t1622 + 0.24000e2_dp &
2138 0 : *t1701*t772*t789 + 0.80e1_dp*t1836 + t1346 + t1838
2139 0 : t1850 = t806*t254
2140 0 : t1853 = t356*t305*t789
2141 0 : t1858 = t802*t142
2142 : t1868 = -0.80e1_dp*t143*t126*t337 + 0.240e2_dp*t755*t371*aa* &
2143 : t757 - 0.2400e2_dp*t1697*t133*t752 - 0.1200e2_dp*t1850 + 0.1200e2_dp &
2144 : *t1853 + 0.80e1_dp*t1657 + t1664 - t1666 - 0.40e1_dp*t1365 + t1570 &
2145 : + t1690 + 0.2e1_dp*t1696 + 0.40e1_dp*t1858*t257 - 0.24000e2_dp*t1119 &
2146 0 : *t133*t90*t91*t789 + 0.80e1_dp*t1671*t327
2147 0 : t1870 = t804*t254
2148 0 : t1872 = t143*t337
2149 0 : t1884 = t1826*t12
2150 0 : t1886 = t1671*t309
2151 0 : t1891 = t808*t254
2152 0 : t1893 = t803*t305
2153 0 : t1894 = t1113*t789
2154 0 : t1897 = t1858*t133
2155 0 : t1901 = t90*t91*t793
2156 : t1904 = -0.80e1_dp*t1870 + 0.80e1_dp*t1872*t92 - 0.160e2_dp*t732*t1 &
2157 : *t758 + 0.1200e2_dp*t771*t772*t793 + 0.2400e2_dp*t1697*af* &
2158 : t773 + t1884*t126 - 0.80e1_dp*t1886*t92 + 0.40e1_dp*t352*t832 &
2159 : - 0.40e1_dp*t1891 + t1893 + t1580 - 0.1200e2_dp*t1894*t92 - 0.40e1_dp &
2160 0 : *t1897*t92 - 0.1200e2_dp*t750*t1901
2161 : t1939 = -t1884 - 0.80e1_dp*t1886 - 0.1200e2_dp*t1894 - 0.40e1_dp*t1829 &
2162 : - 0.40e1_dp*t1897 - 0.2400e2_dp*t1697*t735 - 0.160e2_dp*t732*t338 &
2163 : - 0.24000e2_dp*t1119*t133*t789 - 0.4800e2_dp*t356*t338*t309 &
2164 : - 0.1200e2_dp*t356*t133*t793 + 0.80e1_dp*t1872 + 0.240e2_dp*t143 &
2165 0 : *t1094
2166 : t1945 = -t1730 - 0.240e2_dp*t777*t778*t371 + t1754 - 0.4800e2_dp* &
2167 : t356*t338*t752 + 0.160e2_dp*t732*af*t779 + t1590 + t1764 + &
2168 : 0.2e1_dp*t1765 - t1768 + 0.40e1_dp*t1593 + 0.4800e2_dp*t771*t778* &
2169 : t1804 + t1939*ap*t54 + 0.1200e2_dp*t1113*t828 + t1775 - 0.80e1_dp &
2170 0 : *t1776
2171 : t1949 = t842 - t1609 + t815 + t1611 + t851 + t1782 + t1783 + t597 + &
2172 : t857 + (t1839 + t1868 + t1904 + t1945)*t2 - t817 + t1785 + t1786 &
2173 0 : - t820 + t486 + t603 - t822
2174 :
2175 0 : e_abb(ip) = e_abb(ip) + (t1795 + t1949)*sc
2176 :
2177 0 : t1975 = t1697*t789
2178 0 : t1979 = t789*t309
2179 0 : t1986 = -0.6e1_dp*t337 - 0.6e1_dp*t1094
2180 : t1998 = -0.5699736440e0_dp*t1088*t1979 + 0.2564881396e1_dp*t332*t309 &
2181 : *t793 + 0.2564881399e1_dp*t129*t1986 - 0.5699736440e0_dp*t1099 &
2182 : *t796*t312 + 0.2564881396e1_dp*t343*t312*t799 - 0.2564881399e1_dp &
2183 0 : *t136*t1986
2184 0 : t1999 = t1998*t12
2185 : t2010 = t1085 + 0.120e2_dp*t1797 + 0.3e1_dp*t1613 - 0.3600e2_dp*t356* &
2186 : t309*t1901 + 0.3600e2_dp*t771*t125*t309*t793 + 0.3e1_dp*t1832 &
2187 : + 0.240e2_dp*t1836 - 0.3600e2_dp*t1975*t92 + t1346 + 0.3e1_dp*t1838 &
2188 : + t1999*t126 + 0.40e1_dp*t143*t126*t1986 - 0.3600e2_dp*t1850 &
2189 : + 0.3600e2_dp*t1853 + 0.120e2_dp*t1657 + 0.24000e2_dp*t1119*t126 &
2190 0 : *t1979
2191 0 : t2011 = t1858*t309
2192 0 : t2014 = t732*t793
2193 0 : t2019 = t143*t1986
2194 0 : t2025 = t1119*t1979
2195 : t2048 = -0.120e2_dp*t2011*t92 + t1570 - 0.120e2_dp*t2014*t92 + 0.3e1_dp &
2196 : *t1696 - 0.240e2_dp*t1870 - 0.40e1_dp*t2019*t92 + (-t1999 - 0.120e2_dp &
2197 : *t2011 - 0.3600e2_dp*t1975 - 0.120e2_dp*t2014 - 0.24000e2_dp* &
2198 : t2025 - 0.3600e2_dp*t356*t309*t793 - 0.40e1_dp*t2019)*ap*t54 &
2199 : - 0.24000e2_dp*t2025*t92 - 0.120e2_dp*t1891 + 0.3e1_dp*t1893 + 0.3600e2_dp &
2200 : *t1697*t828 + 0.120e2_dp*t732*t832 + t1998*t56*t92 + &
2201 0 : 0.3e1_dp*t1765 + 0.120e2_dp*t1858*t327 - 0.120e2_dp*t1776
2202 :
2203 : e_bbb(ip) = e_bbb(ip) + (0.3e1_dp*t812 + 0.6e1_dp*t743 + t842 + t851 + t857 + 0.6e1_dp*t762 &
2204 : - 0.240e2_dp*t764 + 0.6e1_dp*t782 + 0.240e2_dp*t784 + 0.3e1_dp*t815 &
2205 : - 0.240e2_dp*t816 - 0.3600e2_dp*t819 - 0.120e2_dp*t821 + 0.3e1_dp*t824 &
2206 : + 0.240e2_dp*t825 + 0.3600e2_dp*t829 + 0.120e2_dp*t833 + (t2010 + &
2207 0 : t2048)*t2)*sc
2208 :
2209 : END IF
2210 : END IF
2211 : END DO
2212 :
2213 : !$OMP END DO
2214 :
2215 2 : END SUBROUTINE vwn5_lsd_calc
2216 :
2217 : END MODULE xc_vwn
2218 :
|
