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