(git:374b731)
Loading...
Searching...
No Matches
xc_vwn.F
Go to the documentation of this file.
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! **************************************************************************************************
21MODULE 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
59CONTAINS
60
61! **************************************************************************************************
62!> \brief ...
63!> \param cutoff ...
64!> \param vwn_params ...
65! **************************************************************************************************
66 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 CALL section_vals_val_get(vwn_params, "functional_type", i_val=type)
74 eps_rho = cutoff
75 CALL set_util(cutoff)
76
77 CALL cite_reference(vosko1980)
78
79 SELECT CASE (type)
80 CASE (do_vwn5)
81 b = 3.72744_dp
82 c = 12.9352_dp
83 x0 = -0.10498_dp
84 bf = 7.06042_dp
85 cf = 18.0578_dp
86 x0f = -0.32500_dp
87 ba = 1.13107_dp
88 ca = 13.0045_dp
89 x0a = -0.00475840_dp
90 CASE (do_vwn3)
91 b = 13.0720_dp
92 c = 42.7198_dp
93 x0 = -0.409286_dp
94 bf = 20.1231_dp
95 cf = 101.578_dp
96 x0f = -0.743294_dp
97 ba = 1.13107_dp
98 ca = 13.0045_dp
99 x0a = -0.00475840_dp
100 CASE DEFAULT
101 cpabort(" Only functionals VWN3 and VWN5 are supported")
102
103 END SELECT
104
105 END SUBROUTINE vwn_init
106
107! **************************************************************************************************
108!> \brief ...
109!> \param reference ...
110!> \param shortform ...
111!> \param needs ...
112!> \param max_deriv ...
113! **************************************************************************************************
114 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 IF (PRESENT(reference)) THEN
120 reference = "S. H. Vosko, L. Wilk and M. Nusair,"// &
121 " Can. J. Phys. 58, 1200 (1980) {LDA version}"
122 END IF
123 IF (PRESENT(shortform)) THEN
124 shortform = "Vosko-Wilk-Nusair Functional {LDA}"
125 END IF
126 IF (PRESENT(needs)) THEN
127 needs%rho = .true.
128 END IF
129 IF (PRESENT(max_deriv)) max_deriv = 3
130
131 END SUBROUTINE vwn_lda_info
132
133! **************************************************************************************************
134!> \brief ...
135!> \param reference ...
136!> \param shortform ...
137!> \param needs ...
138!> \param max_deriv ...
139! **************************************************************************************************
140 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 IF (PRESENT(reference)) THEN
146 reference = "S. H. Vosko, L. Wilk and M. Nusair,"// &
147 " Can. J. Phys. 58, 1200 (1980) {LDA version}"
148 END IF
149 IF (PRESENT(shortform)) THEN
150 shortform = "Vosko-Wilk-Nusair Functional {LDA}"
151 END IF
152 IF (PRESENT(needs)) THEN
153 needs%rho_spin = .true.
154 END IF
155 IF (PRESENT(max_deriv)) max_deriv = 3
156
157 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: x
178 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
179 POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, rho
180 TYPE(xc_derivative_type), POINTER :: deriv
181
182 CALL timeset(routinen, handle)
183
184 CALL section_vals_val_get(vwn_params, "scale_c", r_val=sc)
185
186 CALL xc_rho_set_get(rho_set, rho=rho, &
187 local_bounds=bo, rho_cutoff=epsilon_rho)
188 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
189 CALL vwn_init(epsilon_rho, vwn_params)
190
191 ALLOCATE (x(npoints))
192 CALL calc_srs_pw(rho, x, npoints)
193
194 IF (order >= 1) THEN
195 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
196 allocate_deriv=.true.)
197 CALL xc_derivative_get(deriv, deriv_data=e_0)
198 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
199 allocate_deriv=.true.)
200 CALL xc_derivative_get(deriv, deriv_data=e_rho)
201
202 CALL vwn_lda_01(rho, x, e_0, e_rho, npoints, sc)
203 ELSEIF (order >= 0) THEN
204 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
205 allocate_deriv=.true.)
206 CALL xc_derivative_get(deriv, deriv_data=e_0)
207
208 CALL vwn_lda_0(rho, x, e_0, npoints, sc)
209 ELSEIF (order == -1) THEN
210 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
211 allocate_deriv=.true.)
212 CALL xc_derivative_get(deriv, deriv_data=e_rho)
213
214 CALL vwn_lda_1(rho, x, e_rho, npoints, sc)
215 END IF
216 IF (order >= 2 .OR. order == -2) THEN
217 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
218 allocate_deriv=.true.)
219 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
220
221 CALL vwn_lda_2(rho, x, e_rho_rho, npoints, sc)
222 END IF
223 IF (order >= 3 .OR. order == -3) THEN
224 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
225 allocate_deriv=.true.)
226 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
227
228 CALL vwn_lda_3(rho, x, e_rho_rho_rho, npoints, sc)
229 END IF
230 IF (order > 3 .OR. order < -3) THEN
231 cpabort("derivatives bigger than 3 not implemented")
232 END IF
233
234 DEALLOCATE (x)
235 CALL timestop(handle)
236
237 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 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 q = sqrt(4.0_dp*c - b*b)
258 xb = 2.0_dp*x0 + b
259 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!$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 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 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 q = sqrt(4.0_dp*c - b*b)
303 xb = 2.0_dp*x0 + b
304 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!$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 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 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 q = sqrt(4.0_dp*c - b*b)
355 xb = 2.0_dp*x0 + b
356 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!$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 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 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 q = sqrt(4.0_dp*c - b*b)
408 xb = 2.0_dp*x0 + b
409 px0 = x0*x0 + b*x0 + c
410 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!$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 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 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 q = sqrt(4.0_dp*c - b*b)
468 xb = 2.0_dp*x0 + b
469 px0 = x0*x0 + b*x0 + c
470 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!$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 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 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 POINTER :: dummy, e_0, e_a, e_aa, e_aaa, e_aab, &
537 e_ab, e_abb, e_b, e_bb, e_bbb, rhoa, &
538 rhob
539 TYPE(xc_derivative_type), POINTER :: deriv
540
541 CALL timeset(routinen, handle)
542
543 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 local_bounds=bo, rho_cutoff=epsilon_rho)
547 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
548 CALL vwn_init(epsilon_rho, vwn_params)
549
550 dummy => rhoa
551
552 e_0 => dummy
553 e_a => dummy
554 e_b => dummy
555 e_aa => dummy
556 e_bb => dummy
557 e_ab => dummy
558 e_aaa => dummy
559 e_bbb => dummy
560 e_aab => dummy
561 e_abb => dummy
562
563 IF (order >= 0) THEN
564 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
565 allocate_deriv=.true.)
566 CALL xc_derivative_get(deriv, deriv_data=e_0)
567 END IF
568 IF (order >= 1 .OR. order == -1) THEN
569 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
570 allocate_deriv=.true.)
571 CALL xc_derivative_get(deriv, deriv_data=e_a)
572
573 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
574 allocate_deriv=.true.)
575 CALL xc_derivative_get(deriv, deriv_data=e_b)
576 END IF
577 IF (order >= 2 .OR. order == -2) THEN
578 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
579 allocate_deriv=.true.)
580 CALL xc_derivative_get(deriv, deriv_data=e_aa)
581
582 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
583 allocate_deriv=.true.)
584 CALL xc_derivative_get(deriv, deriv_data=e_bb)
585
586 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
587 allocate_deriv=.true.)
588 CALL xc_derivative_get(deriv, deriv_data=e_ab)
589 END IF
590 IF (order >= 3 .OR. order == -3) THEN
591 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
592 allocate_deriv=.true.)
593 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 allocate_deriv=.true.)
597 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 allocate_deriv=.true.)
601 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 allocate_deriv=.true.)
605 CALL xc_derivative_get(deriv, deriv_data=e_abb)
606
607 END IF
608 IF (order > 3 .OR. order < -3) THEN
609 cpabort("derivatives bigger than 3 not implemented")
610 END IF
611
612 CALL section_vals_val_get(vwn_params, "functional_type", i_val=type)
613 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!$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!$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 cpabort(" Only functionals VWN3 and VWN5 are supported")
650 END SELECT
651
652 CALL timestop(handle)
653
654 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 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 cp = c
709 ap = a
710 bp = b
711 x0p = x0
712 qp = sqrt(4.0_dp*cp - bp*bp)
713 qf = sqrt(4.0_dp*cf - bf*bf)
714
715!$OMP DO
716 DO ip = 1, npoints
717 myrhoa = max(rhoa(ip), 0.0_dp)
718 myrhob = max(rhob(ip), 0.0_dp)
719 myrho = myrhoa + myrhob
720 IF (myrho > eps_rho) THEN
721 myrhoa = max(epsilon(0.0_dp)*1.e4_dp, myrhoa)
722 myrhob = max(epsilon(0.0_dp)*1.e4_dp, myrhob)
723 myrho = myrhoa + myrhob
724
725 IF (order >= 0) THEN
726 t1 = 0.1e1_dp/0.3141592654e1_dp
727 t2 = myrho
728 t3 = 0.1e1_dp/t2
729 t4 = t1*t3
730 t5 = t4**0.3333333332e0_dp
731 t6 = 0.9085602964e0_dp*t5
732 t7 = t4**0.1666666666e0_dp
733 t10 = t6 + 0.9531842930e0_dp*bp*t7 + cp
734 t11 = 0.1e1_dp/t10
735 t14 = log(0.9085602964e0_dp*t5*t11)
736 t15 = 0.1906368586e1_dp*t7
737 t16 = t15 + bp
738 t19 = atan(qp/t16)
739 t21 = 0.1e1_dp/qp
740 t24 = bp*x0p
741 t25 = 0.9531842930e0_dp*t7
742 t26 = t25 - x0p
743 t27 = t26**0.20e1_dp
744 t29 = log(t27*t11)
745 t32 = 0.20e1_dp*bp + 0.40e1_dp*x0p
746 t36 = x0p**0.20e1_dp
747 t38 = 0.1e1_dp/(t36 + t24 + cp)
748 t42 = ap*(t14 + 0.20e1_dp*bp*t19*t21 - t24*(t29 + t32*t19 &
749 *t21)*t38)
750 t43 = myrhoa - myrhob
751 t44 = t43*t3
752 t45 = 0.10e1_dp + t44
753 t46 = t45**0.1333333333e1_dp
754 t48 = 0.10e1_dp - t44
755 t49 = t48**0.1333333333e1_dp
756 t51 = 0.1923661050e1_dp*t46 + 0.1923661050e1_dp*t49 - 0.3847322100e1_dp
757 t54 = t6 + 0.9531842930e0_dp*bf*t7 + cf
758 t55 = 0.1e1_dp/t54
759 t58 = log(0.9085602964e0_dp*t5*t55)
760 t59 = t15 + bf
761 t62 = atan(qf/t59)
762 t64 = 0.1e1_dp/qf
763 t67 = bf*x0f
764 t68 = t25 - x0f
765 t69 = t68**0.20e1_dp
766 t71 = log(t69*t55)
767 t74 = 0.20e1_dp*bf + 0.40e1_dp*x0f
768 t78 = x0f**0.20e1_dp
769 t80 = 0.1e1_dp/(t78 + t67 + cf)
770 t85 = af*(t58 + 0.20e1_dp*bf*t62*t64 - t67*(t71 + t74*t62 &
771 *t64)*t80) - t42
772 t86 = t51*t85
773
774 e_0(ip) = e_0(ip) + ((t42 + t86)*t2)*sc
775
776 END IF
777 IF (order >= 1 .OR. order == -1) THEN
778 t88 = t4**(-0.6666666668e0_dp)
779 t89 = t88*t11
780 t90 = t2**2
781 t91 = 0.1e1_dp/t90
782 t92 = t1*t91
783 t95 = t10**2
784 t96 = 0.1e1_dp/t95
785 t97 = t5*t96
786 t98 = t88*t1
787 t100 = 0.3028534320e0_dp*t98*t91
788 t101 = t4**(-0.8333333334e0_dp)
789 t102 = bp*t101
790 t105 = -t100 - 0.1588640488e0_dp*t102*t92
791 t108 = -0.3028534320e0_dp*t89*t92 - 0.9085602964e0_dp*t97*t105
792 t109 = t4**(-0.3333333332e0_dp)
793 t110 = t108*t109
794 t113 = t16**2
795 t114 = 0.1e1_dp/t113
796 t115 = bp*t114
797 t116 = t115*t101
798 t117 = qp**2
799 t119 = 0.1e1_dp + t117*t114
800 t120 = 0.1e1_dp/t119
801 t121 = t92*t120
802 t124 = t26**0.10e1_dp
803 t125 = t124*t11
804 t126 = t101*t1
805 t127 = t126*t91
806 t130 = t27*t96
807 t132 = -0.3177280976e0_dp*t125*t127 - t130*t105
808 t133 = t26**(-0.20e1_dp)
809 t134 = t132*t133
810 t136 = t32*t114
811 t137 = t136*t101
812 t144 = ap*(0.1100642416e1_dp*t110*t10 + 0.6354561950e0_dp*t116* &
813 t121 - t24*(t134*t10 + 0.3177280975e0_dp*t137*t121)*t38)
814 t145 = t45**0.333333333e0_dp
815 t146 = t43*t91
816 t147 = t3 - t146
817 t150 = t48**0.333333333e0_dp
818 t151 = -t147
819 t154 = 0.2564881399e1_dp*t145*t147 + 0.2564881399e1_dp*t150*t151
820 t155 = t154*t85
821 t156 = t88*t55
822 t159 = t54**2
823 t160 = 0.1e1_dp/t159
824 t161 = t5*t160
825 t162 = bf*t101
826 t165 = -t100 - 0.1588640488e0_dp*t162*t92
827 t168 = -0.3028534320e0_dp*t156*t92 - 0.9085602964e0_dp*t161*t165
828 t169 = t168*t109
829 t172 = t59**2
830 t173 = 0.1e1_dp/t172
831 t174 = bf*t173
832 t175 = t174*t101
833 t176 = qf**2
834 t178 = 0.1e1_dp + t176*t173
835 t179 = 0.1e1_dp/t178
836 t180 = t92*t179
837 t183 = t68**0.10e1_dp
838 t184 = t183*t55
839 t187 = t69*t160
840 t189 = -0.3177280976e0_dp*t184*t127 - t187*t165
841 t190 = t68**(-0.20e1_dp)
842 t191 = t189*t190
843 t193 = t74*t173
844 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 t144
848 t203 = t51*t202
849
850 e_a(ip) = e_a(ip) + ((t144 + t155 + t203)*t2 + t42 + t86)*sc
851
852 t206 = -t3 - t146
853 t209 = -t206
854 t212 = 0.2564881399e1_dp*t145*t206 + 0.2564881399e1_dp*t150*t209
855 t213 = t212*t85
856
857 e_b(ip) = e_b(ip) + ((t144 + t213 + t203)*t2 + t42 + t86)*sc
858
859 END IF
860 IF (order >= 2 .OR. order == -2) THEN
861 t216 = t4**(-0.1666666667e1_dp)
862 t217 = t216*t11
863 t218 = 0.3141592654e1_dp**2
864 t219 = 0.1e1_dp/t218
865 t220 = t90**2
866 t221 = 0.1e1_dp/t220
867 t222 = t219*t221
868 t223 = t217*t222
869 t225 = t88*t96
870 t226 = t92*t105
871 t230 = 0.1e1_dp/t90/t2
872 t231 = t1*t230
873 t235 = 0.1e1_dp/t95/t10
874 t236 = t5*t235
875 t237 = t105**2
876 t240 = t216*t219
877 t241 = t240*t221
878 t242 = 0.2019022880e0_dp*t241
879 t244 = 0.6057068640e0_dp*t98*t230
880 t245 = t4**(-0.1833333333e1_dp)
881 t246 = bp*t245
882 t251 = -t242 + t244 - 0.1323867073e0_dp*t246*t222 + 0.3177280976e0_dp &
883 *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 *t97*t251
887 t255 = t254*t109
888 t258 = t4**(-0.1333333333e1_dp)
889 t259 = t108*t258
890 t260 = t10*t1
891 t261 = t260*t91
892 t267 = 0.1e1_dp/t113/t16
893 t268 = bp*t267
894 t269 = t268*t216
895 t270 = t222*t120
896 t273 = t115*t245
897 t276 = t231*t120
898 t279 = t113**2
899 t281 = 0.1e1_dp/t279/t16
900 t282 = bp*t281
901 t283 = t282*t216
902 t284 = t119**2
903 t286 = 0.1e1_dp/t284*t117
904 t287 = t222*t286
905 t291 = t124*t96
906 t292 = t291*t101
907 t295 = t245*t219
908 t296 = t295*t221
909 t299 = t126*t230
910 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 t302*t237 - t130*t251
914 t307 = t306*t133
915 t309 = t26**(-0.30e1_dp)
916 t310 = t132*t309
917 t311 = t310*t10
918 t315 = t32*t267
919 t316 = t315*t216
920 t319 = t136*t245
921 t324 = t32*t281
922 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 *t325*t287)*t38)
930 t333 = t45**(-0.666666667e0_dp)
931 t334 = t147**2
932 t337 = t43*t230
933 t339 = -0.2e1_dp*t91 + 0.2e1_dp*t337
934 t342 = t48**(-0.666666667e0_dp)
935 t343 = t151**2
936 t346 = -t339
937 t349 = 0.8549604655e0_dp*t333*t334 + 0.2564881399e1_dp*t145*t339 &
938 + 0.8549604655e0_dp*t342*t343 + 0.2564881399e1_dp*t150*t346
939 t350 = t349*t85
940 t351 = t154*t202
941 t352 = 0.2e1_dp*t351
942 t353 = t216*t55
943 t354 = t353*t222
944 t356 = t88*t160
945 t357 = t92*t165
946 t363 = 0.1e1_dp/t159/t54
947 t364 = t5*t363
948 t365 = t165**2
949 t368 = bf*t245
950 t373 = -t242 + t244 - 0.1323867073e0_dp*t368*t222 + 0.3177280976e0_dp &
951 *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 *t161*t373
955 t377 = t376*t109
956 t380 = t168*t258
957 t381 = t54*t1
958 t382 = t381*t91
959 t388 = 0.1e1_dp/t172/t59
960 t389 = bf*t388
961 t390 = t389*t216
962 t391 = t222*t179
963 t394 = t174*t245
964 t397 = t231*t179
965 t400 = t172**2
966 t402 = 0.1e1_dp/t400/t59
967 t403 = bf*t402
968 t404 = t403*t216
969 t405 = t178**2
970 t407 = 0.1e1_dp/t405*t176
971 t408 = t222*t407
972 t412 = t183*t160
973 t413 = t412*t101
974 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 t420*t365 - t187*t373
978 t425 = t424*t190
979 t427 = t68**(-0.30e1_dp)
980 t428 = t189*t427
981 t429 = t428*t54
982 t433 = t74*t388
983 t434 = t433*t216
984 t437 = t193*t245
985 t442 = t74*t402
986 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 *t443*t408)*t80) - t332
994 t452 = t51*t451
995 t455 = 0.2e1_dp*t144
996 t457 = 0.2e1_dp*t203
997
998 e_aa(ip) = e_aa(ip) + ((t332 + t350 + t352 + t452)*t2 + t455 + 0.2e1_dp*t155 + t457)*sc
999
1000 t458 = t333*t147
1001 t461 = t145*t43
1002 t464 = t342*t151
1003 t467 = t150*t43
1004 t470 = 0.8549604655e0_dp*t458*t206 + 0.5129762798e1_dp*t461*t230 &
1005 + 0.8549604655e0_dp*t464*t209 - 0.5129762798e1_dp*t467*t230
1006 t471 = t470*t85
1007 t472 = t212*t202
1008
1009 e_ab(ip) = e_ab(ip) + ((t332 + t471 + t351 + t472 + t452)*t2 + t455 + t155 + t457 &
1010 + t213)*sc
1011 t475 = t206**2
1012 t479 = 0.2e1_dp*t91 + 0.2e1_dp*t337
1013 t482 = t209**2
1014 t485 = -t479
1015 t488 = 0.8549604655e0_dp*t333*t475 + 0.2564881399e1_dp*t145*t479 &
1016 + 0.8549604655e0_dp*t342*t482 + 0.2564881399e1_dp*t150*t485
1017 t489 = t488*t85
1018 t490 = 0.2e1_dp*t472
1019
1020 e_bb(ip) = e_bb(ip) + ((t332 + t489 + t490 + t452)*t2 + t455 + 0.2e1_dp*t213 + t457)*sc
1021
1022 END IF
1023 IF (order >= 3 .OR. order == -3) THEN
1024 t494 = t4**(-0.2666666667e1_dp)
1025 t497 = 0.1e1_dp/t218/0.3141592654e1_dp
1026 t499 = 0.1e1_dp/t220/t90
1027 t500 = t497*t499
1028 t501 = t494*t11*t500
1029 t504 = t222*t105
1030 t505 = t216*t96*t504
1031 t508 = 0.1e1_dp/t220/t2
1032 t509 = t219*t508
1033 t510 = t217*t509
1034 t513 = t92*t237
1035 t516 = t231*t105
1036 t519 = t92*t251
1037 t522 = t1*t221
1038 t525 = t95**2
1039 t526 = 0.1e1_dp/t525
1040 t528 = t237*t105
1041 t531 = t105*t251
1042 t536 = 0.3365038134e0_dp*t494*t497*t499
1043 t538 = 0.1211413728e1_dp*t240*t508
1044 t540 = 0.1817120592e1_dp*t98*t221
1045 t541 = t4**(-0.2833333333e1_dp)
1046 t549 = -t536 + t538 - t540 - 0.2427089633e0_dp*bp*t541*t500 + 0.7943202439e0_dp &
1047 *t246*t509 - 0.9531842928e0_dp*t102*t522
1048 t559 = t500*t120
1049 t563 = 0.1e1_dp/t279/t113
1050 t565 = t4**(-0.2500000000e1_dp)
1051 t567 = t500*t286
1052 t570 = t509*t120
1053 t573 = t4**(-0.2666666666e1_dp)
1054 t583 = t522*t120
1055 t589 = t509*t286
1056 t595 = t279**2
1057 t596 = 0.1e1_dp/t595
1058 t601 = t117**2
1059 t603 = t500/t284/t119*t601
1060 t613 = t4**(-0.2333333333e1_dp)
1061 t625 = 0.1e1_dp/t279
1062 t665 = t26**(-0.40e1_dp)
1063 t690 = t541*t497*t499
1064 t693 = t295*t508
1065 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 *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 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 *bp*t625*t565*t559
1100 t720 = ap*t719
1101 t721 = t45**(-0.1666666667e1_dp)
1102 t727 = t43*t221
1103 t729 = 0.6e1_dp*t230 - 0.6e1_dp*t727
1104 t732 = t48**(-0.1666666667e1_dp)
1105 t743 = t349*t202
1106 t745 = t154*t451
1107 t748 = t500*t179
1108 t752 = t500*t407
1109 t755 = t522*t179
1110 t758 = t509*t407
1111 t769 = t68**(-0.40e1_dp)
1112 t792 = t400**2
1113 t793 = 0.1e1_dp/t792
1114 t798 = t176**2
1115 t800 = t500/t405/t178*t798
1116 t804 = t494*t55*t500
1117 t807 = t222*t165
1118 t808 = t216*t160*t807
1119 t810 = t353*t509
1120 t814 = t92*t365
1121 t820 = t231*t165
1122 t823 = t92*t373
1123 t835 = t159**2
1124 t836 = 0.1e1_dp/t835
1125 t838 = t365*t165
1126 t841 = t165*t373
1127 t851 = -t536 + t538 - t540 - 0.2427089633e0_dp*bf*t541*t500 + 0.7943202439e0_dp &
1128 *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 *t841 - t187*t851
1136 t860 = t509*t179
1137 t863 = 0.1e1_dp/t400
1138 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.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 *t389*t573*t748
1166 t950 = t51*(af*t947 - t720)
1167 t953 = 0.3e1_dp*t332
1168 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 + t953 + 0.3e1_dp*t350 + 0.6e1_dp*t351 + t956)*sc
1175
1176 t967 = 0.2e1_dp*t230 - 0.6e1_dp*t727
1177 t984 = 0.2e1_dp*t470*t202
1178 t986 = t212*t451
1179 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 + t950)*t2 + t953 + t350 + 0.4e1_dp*t351 + t956 + t990 + t490)*sc
1187
1188 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 + t489)*sc
1198
1199 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 + 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 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 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 cp = c
1296 ap = a
1297 bp = b
1298 x0p = x0
1299 qp = sqrt(4.0_dp*cp - bp*bp)
1300 qa = sqrt(4.0_dp*ca - ba*ba)
1301 qf = sqrt(4.0_dp*cf - bf*bf)
1302 d2f0 = 4.0_dp/(9.0_dp*(2.0_dp**(1.0_dp/3.0_dp) - 1.0_dp))
1303
1304!$OMP DO
1305
1306 DO ip = 1, npoints
1307 myrhoa = max(rhoa(ip), 0.0_dp)
1308 myrhob = max(rhob(ip), 0.0_dp)
1309 myrho = myrhoa + myrhob
1310 IF (myrho > eps_rho) THEN
1311 myrhoa = max(epsilon(0.0_dp)*1.e4_dp, myrhoa)
1312 myrhob = max(epsilon(0.0_dp)*1.e4_dp, myrhob)
1313 myrho = myrhoa + myrhob
1314
1315 IF (order >= 0) THEN
1316 t1 = myrhoa - myrhob
1317 t2 = myrho
1318 t3 = 0.1e1_dp/t2
1319 t4 = t1*t3
1320 t5 = 0.10e1_dp + t4
1321 t6 = t5**0.1333333333e1_dp
1322 t8 = 0.10e1_dp - t4
1323 t9 = t8**0.1333333333e1_dp
1324 t11 = 0.1923661050e1_dp*t6 + 0.1923661050e1_dp*t9 - 0.3847322100e1_dp
1325 t12 = t4**0.40e1_dp
1326 t13 = t11*t12
1327 t15 = (0.10e1_dp - t13)*ap
1328 t16 = 0.1e1_dp/0.3141592654e1_dp
1329 t17 = t16*t3
1330 t18 = t17**0.3333333332e0_dp
1331 t19 = 0.9085602964e0_dp*t18
1332 t20 = t17**0.1666666666e0_dp
1333 t23 = t19 + 0.9531842930e0_dp*bp*t20 + cp
1334 t24 = 0.1e1_dp/t23
1335 t27 = log(0.9085602964e0_dp*t18*t24)
1336 t28 = 0.1906368586e1_dp*t20
1337 t29 = t28 + bp
1338 t32 = atan(qp/t29)
1339 t34 = 0.1e1_dp/qp
1340 t37 = bp*x0p
1341 t38 = 0.9531842930e0_dp*t20
1342 t39 = t38 - x0p
1343 t40 = t39**0.20e1_dp
1344 t42 = log(t40*t24)
1345 t45 = 0.20e1_dp*bp + 0.40e1_dp*x0p
1346 t49 = x0p**0.20e1_dp
1347 t51 = 0.1e1_dp/(t49 + t37 + cp)
1348 t54 = t27 + 0.20e1_dp*bp*t32*t34 - t37*(t42 + t45*t32*t34) &
1349 *t51
1350 t55 = t15*t54
1351 t56 = 0.10e1_dp - t12
1352 t57 = t11*t56
1353 t60 = t19 + 0.9531842930e0_dp*ba*t20 + ca
1354 t61 = 0.1e1_dp/t60
1355 t64 = log(0.9085602964e0_dp*t18*t61)
1356 t65 = t28 + ba
1357 t68 = atan(qa/t65)
1358 t70 = 0.1e1_dp/qa
1359 t73 = ba*x0a
1360 t74 = t38 - x0a
1361 t75 = t74**0.20e1_dp
1362 t77 = log(t75*t61)
1363 t80 = 0.20e1_dp*ba + 0.40e1_dp*x0a
1364 t84 = x0a**0.20e1_dp
1365 t86 = 0.1e1_dp/(t84 + t73 + ca)
1366 t89 = t64 + 0.20e1_dp*ba*t68*t70 - t73*(t77 + t80*t68*t70) &
1367 *t86
1368 t90 = aa*t89
1369 t91 = 0.1e1_dp/d2f0
1370 t92 = t90*t91
1371 t93 = t57*t92
1372 t96 = t19 + 0.9531842930e0_dp*bf*t20 + cf
1373 t97 = 0.1e1_dp/t96
1374 t100 = log(0.9085602964e0_dp*t18*t97)
1375 t101 = t28 + bf
1376 t104 = atan(qf/t101)
1377 t106 = 0.1e1_dp/qf
1378 t109 = bf*x0f
1379 t110 = t38 - x0f
1380 t111 = t110**0.20e1_dp
1381 t113 = log(t111*t97)
1382 t116 = 0.20e1_dp*bf + 0.40e1_dp*x0f
1383 t120 = x0f**0.20e1_dp
1384 t122 = 0.1e1_dp/(t120 + t109 + cf)
1385 t125 = t100 + 0.20e1_dp*bf*t104*t106 - t109*(t113 + t116*t104 &
1386 *t106)*t122
1387 t126 = af*t125
1388 t127 = t13*t126
1389
1390 e_0(ip) = e_0(ip) + ((t55 + t93 + t127)*t2)*sc
1391
1392 END IF
1393 IF (order >= 1 .OR. order == -1) THEN
1394 t129 = t5**0.333333333e0_dp
1395 t130 = t2**2
1396 t131 = 0.1e1_dp/t130
1397 t132 = t1*t131
1398 t133 = t3 - t132
1399 t136 = t8**0.333333333e0_dp
1400 t137 = -t133
1401 t140 = 0.2564881399e1_dp*t129*t133 + 0.2564881399e1_dp*t136*t137
1402 t141 = t140*t12
1403 t142 = t4**0.30e1_dp
1404 t143 = t11*t142
1405 t144 = t143*t133
1406 t147 = (-t141 - 0.40e1_dp*t144)*ap
1407 t148 = t147*t54
1408 t149 = t17**(-0.6666666668e0_dp)
1409 t150 = t149*t24
1410 t151 = t16*t131
1411 t154 = t23**2
1412 t155 = 0.1e1_dp/t154
1413 t156 = t18*t155
1414 t157 = t149*t16
1415 t159 = 0.3028534320e0_dp*t157*t131
1416 t160 = t17**(-0.8333333334e0_dp)
1417 t161 = bp*t160
1418 t164 = -t159 - 0.1588640488e0_dp*t161*t151
1419 t167 = -0.3028534320e0_dp*t150*t151 - 0.9085602964e0_dp*t156*t164
1420 t168 = t17**(-0.3333333332e0_dp)
1421 t169 = t167*t168
1422 t172 = t29**2
1423 t173 = 0.1e1_dp/t172
1424 t174 = bp*t173
1425 t175 = t174*t160
1426 t176 = qp**2
1427 t178 = 0.1e1_dp + t176*t173
1428 t179 = 0.1e1_dp/t178
1429 t180 = t151*t179
1430 t183 = t39**0.10e1_dp
1431 t184 = t183*t24
1432 t185 = t160*t16
1433 t186 = t185*t131
1434 t189 = t40*t155
1435 t191 = -0.3177280976e0_dp*t184*t186 - t189*t164
1436 t192 = t39**(-0.20e1_dp)
1437 t193 = t191*t192
1438 t195 = t45*t173
1439 t196 = t195*t160
1440 t202 = 0.1100642416e1_dp*t169*t23 + 0.6354561950e0_dp*t175*t180 - &
1441 t37*(t193*t23 + 0.3177280975e0_dp*t196*t180)*t51
1442 t203 = t15*t202
1443 t204 = t140*t56
1444 t205 = t204*t92
1445 t206 = t144*t92
1446 t207 = 0.40e1_dp*t206
1447 t208 = t149*t61
1448 t211 = t60**2
1449 t212 = 0.1e1_dp/t211
1450 t213 = t18*t212
1451 t214 = ba*t160
1452 t217 = -t159 - 0.1588640488e0_dp*t214*t151
1453 t220 = -0.3028534320e0_dp*t208*t151 - 0.9085602964e0_dp*t213*t217
1454 t221 = t220*t168
1455 t224 = t65**2
1456 t225 = 0.1e1_dp/t224
1457 t226 = ba*t225
1458 t227 = t226*t160
1459 t228 = qa**2
1460 t230 = 0.1e1_dp + t228*t225
1461 t231 = 0.1e1_dp/t230
1462 t232 = t151*t231
1463 t235 = t74**0.10e1_dp
1464 t236 = t235*t61
1465 t239 = t75*t212
1466 t241 = -0.3177280976e0_dp*t236*t186 - t239*t217
1467 t242 = t74**(-0.20e1_dp)
1468 t243 = t241*t242
1469 t245 = t80*t225
1470 t246 = t245*t160
1471 t252 = 0.1100642416e1_dp*t221*t60 + 0.6354561950e0_dp*t227*t232 - &
1472 t73*(t243*t60 + 0.3177280975e0_dp*t246*t232)*t86
1473 t253 = aa*t252
1474 t254 = t253*t91
1475 t255 = t57*t254
1476 t256 = t141*t126
1477 t257 = t126*t133
1478 t258 = t143*t257
1479 t259 = 0.40e1_dp*t258
1480 t260 = t149*t97
1481 t263 = t96**2
1482 t264 = 0.1e1_dp/t263
1483 t265 = t18*t264
1484 t266 = bf*t160
1485 t269 = -t159 - 0.1588640488e0_dp*t266*t151
1486 t272 = -0.3028534320e0_dp*t260*t151 - 0.9085602964e0_dp*t265*t269
1487 t273 = t272*t168
1488 t276 = t101**2
1489 t277 = 0.1e1_dp/t276
1490 t278 = bf*t277
1491 t279 = t278*t160
1492 t280 = qf**2
1493 t282 = 0.1e1_dp + t280*t277
1494 t283 = 0.1e1_dp/t282
1495 t284 = t151*t283
1496 t287 = t110**0.10e1_dp
1497 t288 = t287*t97
1498 t291 = t111*t264
1499 t293 = -0.3177280976e0_dp*t288*t186 - t291*t269
1500 t294 = t110**(-0.20e1_dp)
1501 t295 = t293*t294
1502 t297 = t116*t277
1503 t298 = t297*t160
1504 t304 = 0.1100642416e1_dp*t273*t96 + 0.6354561950e0_dp*t279*t284 - &
1505 t109*(t295*t96 + 0.3177280975e0_dp*t298*t284)*t122
1506 t305 = af*t304
1507 t306 = t13*t305
1508
1509 e_a(ip) = e_a(ip) + ((t148 + t203 + t205 - t207 + t255 + t256 + t259 + t306)*t2 + &
1510 t55 + t93 + t127)*sc
1511
1512 t309 = -t3 - t132
1513 t312 = -t309
1514 t315 = 0.2564881399e1_dp*t129*t309 + 0.2564881399e1_dp*t136*t312
1515 t316 = t315*t12
1516 t317 = t143*t309
1517 t320 = (-t316 - 0.40e1_dp*t317)*ap
1518 t321 = t320*t54
1519 t322 = t315*t56
1520 t323 = t322*t92
1521 t324 = t317*t92
1522 t325 = 0.40e1_dp*t324
1523 t326 = t316*t126
1524 t327 = t126*t309
1525 t328 = t143*t327
1526 t329 = 0.40e1_dp*t328
1527
1528 e_b(ip) = e_b(ip) + ((t321 + t203 + t323 - t325 + t255 + t326 + t329 + t306)*t2 + &
1529 t55 + t93 + t127)*sc
1530
1531 END IF
1532 IF (order >= 2 .OR. order == -2) THEN
1533 t332 = t5**(-0.666666667e0_dp)
1534 t333 = t133**2
1535 t337 = 0.1e1_dp/t130/t2
1536 t338 = t1*t337
1537 t340 = -0.2e1_dp*t131 + 0.2e1_dp*t338
1538 t343 = t8**(-0.666666667e0_dp)
1539 t344 = t137**2
1540 t347 = -t340
1541 t350 = 0.8549604655e0_dp*t332*t333 + 0.2564881399e1_dp*t129*t340 &
1542 + 0.8549604655e0_dp*t343*t344 + 0.2564881399e1_dp*t136*t347
1543 t351 = t350*t12
1544 t352 = t140*t142
1545 t353 = t352*t133
1546 t355 = t4**0.20e1_dp
1547 t356 = t11*t355
1548 t357 = t356*t333
1549 t359 = t143*t340
1550 t362 = (-t351 - 0.80e1_dp*t353 - 0.1200e2_dp*t357 - 0.40e1_dp*t359)* &
1551 ap
1552 t363 = t362*t54
1553 t364 = t147*t202
1554 t365 = 0.2e1_dp*t364
1555 t366 = t17**(-0.1666666667e1_dp)
1556 t367 = t366*t24
1557 t368 = 0.3141592654e1_dp**2
1558 t369 = 0.1e1_dp/t368
1559 t370 = t130**2
1560 t371 = 0.1e1_dp/t370
1561 t372 = t369*t371
1562 t373 = t367*t372
1563 t375 = t149*t155
1564 t376 = t151*t164
1565 t379 = t16*t337
1566 t383 = 0.1e1_dp/t154/t23
1567 t384 = t18*t383
1568 t385 = t164**2
1569 t388 = t366*t369
1570 t389 = t388*t371
1571 t390 = 0.2019022880e0_dp*t389
1572 t392 = 0.6057068640e0_dp*t157*t337
1573 t393 = t17**(-0.1833333333e1_dp)
1574 t394 = bp*t393
1575 t399 = -t390 + t392 - 0.1323867073e0_dp*t394*t372 + 0.3177280976e0_dp &
1576 *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 *t156*t399
1580 t403 = t402*t168
1581 t406 = t17**(-0.1333333333e1_dp)
1582 t407 = t167*t406
1583 t408 = t23*t16
1584 t409 = t408*t131
1585 t415 = 0.1e1_dp/t172/t29
1586 t416 = bp*t415
1587 t417 = t416*t366
1588 t418 = t372*t179
1589 t421 = t174*t393
1590 t424 = t379*t179
1591 t427 = t172**2
1592 t429 = 0.1e1_dp/t427/t29
1593 t430 = bp*t429
1594 t431 = t430*t366
1595 t432 = t178**2
1596 t434 = 0.1e1_dp/t432*t176
1597 t435 = t372*t434
1598 t439 = t183*t155
1599 t440 = t439*t160
1600 t443 = t393*t369
1601 t444 = t443*t371
1602 t447 = t185*t337
1603 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 t450*t385 - t189*t399
1607 t455 = t454*t192
1608 t457 = t39**(-0.30e1_dp)
1609 t458 = t191*t457
1610 t459 = t458*t23
1611 t463 = t45*t415
1612 t464 = t463*t366
1613 t467 = t195*t393
1614 t472 = t45*t429
1615 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 t473*t435)*t51
1623 t480 = t15*t479
1624 t481 = t350*t56
1625 t482 = t481*t92
1626 t483 = t353*t92
1627 t484 = 0.80e1_dp*t483
1628 t485 = t204*t254
1629 t486 = 0.2e1_dp*t485
1630 t487 = t357*t92
1631 t488 = 0.1200e2_dp*t487
1632 t489 = t359*t92
1633 t490 = 0.40e1_dp*t489
1634 t491 = t144*t254
1635 t492 = 0.80e1_dp*t491
1636 t493 = t366*t61
1637 t494 = t493*t372
1638 t496 = t149*t212
1639 t497 = t151*t217
1640 t503 = 0.1e1_dp/t211/t60
1641 t504 = t18*t503
1642 t505 = t217**2
1643 t508 = ba*t393
1644 t513 = -t390 + t392 - 0.1323867073e0_dp*t508*t372 + 0.3177280976e0_dp &
1645 *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 *t213*t513
1649 t517 = t516*t168
1650 t520 = t220*t406
1651 t521 = t60*t16
1652 t522 = t521*t131
1653 t528 = 0.1e1_dp/t224/t65
1654 t529 = ba*t528
1655 t530 = t529*t366
1656 t531 = t372*t231
1657 t534 = t226*t393
1658 t537 = t379*t231
1659 t540 = t224**2
1660 t542 = 0.1e1_dp/t540/t65
1661 t543 = ba*t542
1662 t544 = t543*t366
1663 t545 = t230**2
1664 t547 = 0.1e1_dp/t545*t228
1665 t548 = t372*t547
1666 t552 = t235*t212
1667 t553 = t552*t160
1668 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 t560*t505 - t239*t513
1672 t565 = t564*t242
1673 t567 = t74**(-0.30e1_dp)
1674 t568 = t241*t567
1675 t569 = t568*t60
1676 t573 = t80*t528
1677 t574 = t573*t366
1678 t577 = t245*t393
1679 t582 = t80*t542
1680 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 *t583*t548)*t86)*t91
1688 t592 = t57*t591
1689 t593 = t351*t126
1690 t594 = t352*t257
1691 t595 = 0.80e1_dp*t594
1692 t596 = t141*t305
1693 t597 = 0.2e1_dp*t596
1694 t598 = t126*t333
1695 t599 = t356*t598
1696 t600 = 0.1200e2_dp*t599
1697 t601 = t305*t133
1698 t602 = t143*t601
1699 t603 = 0.80e1_dp*t602
1700 t604 = t126*t340
1701 t605 = t143*t604
1702 t606 = 0.40e1_dp*t605
1703 t607 = t366*t97
1704 t608 = t607*t372
1705 t610 = t149*t264
1706 t611 = t151*t269
1707 t617 = 0.1e1_dp/t263/t96
1708 t618 = t18*t617
1709 t619 = t269**2
1710 t622 = bf*t393
1711 t627 = -t390 + t392 - 0.1323867073e0_dp*t622*t372 + 0.3177280976e0_dp &
1712 *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 *t265*t627
1716 t631 = t630*t168
1717 t634 = t272*t406
1718 t635 = t96*t16
1719 t636 = t635*t131
1720 t642 = 0.1e1_dp/t276/t101
1721 t643 = bf*t642
1722 t644 = t643*t366
1723 t645 = t372*t283
1724 t648 = t278*t393
1725 t651 = t379*t283
1726 t654 = t276**2
1727 t656 = 0.1e1_dp/t654/t101
1728 t657 = bf*t656
1729 t658 = t657*t366
1730 t659 = t282**2
1731 t661 = 0.1e1_dp/t659*t280
1732 t662 = t372*t661
1733 t666 = t287*t264
1734 t667 = t666*t160
1735 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 t674*t619 - t291*t627
1739 t679 = t678*t294
1740 t681 = t110**(-0.30e1_dp)
1741 t682 = t293*t681
1742 t683 = t682*t96
1743 t687 = t116*t642
1744 t688 = t687*t366
1745 t691 = t297*t393
1746 t696 = t116*t656
1747 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 *t697*t662)*t122)
1755 t705 = t13*t704
1756 t706 = t363 + t365 + t480 + t482 - t484 + t486 - t488 - t490 - t492 &
1757 + t592 + t593 + t595 + t597 + t600 + t603 + t606 + t705
1758 t709 = 0.2e1_dp*t203
1759 t712 = 0.2e1_dp*t255
1760 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 + t712 + 0.2e1_dp*t256 + 0.80e1_dp*t258 + t715)*sc
1764
1765 t716 = t332*t133
1766 t719 = t129*t1
1767 t722 = t343*t137
1768 t725 = t136*t1
1769 t728 = 0.8549604655e0_dp*t716*t309 + 0.5129762798e1_dp*t719*t337 &
1770 + 0.8549604655e0_dp*t722*t312 - 0.5129762798e1_dp*t725*t337
1771 t729 = t728*t12
1772 t730 = t352*t309
1773 t732 = t315*t142
1774 t733 = t732*t133
1775 t735 = t133*t309
1776 t741 = (-t729 - 0.40e1_dp*t730 - 0.40e1_dp*t733 - 0.1200e2_dp*t356*t735 &
1777 - 0.80e1_dp*t143*t338)*ap
1778 t742 = t741*t54
1779 t743 = t320*t202
1780 t744 = t728*t56
1781 t745 = t744*t92
1782 t746 = t730*t92
1783 t748 = t733*t92
1784 t750 = t356*t133
1785 t751 = t91*t309
1786 t752 = t90*t751
1787 t753 = t750*t752
1788 t755 = t143*t1
1789 t756 = t337*aa
1790 t757 = t89*t91
1791 t758 = t756*t757
1792 t759 = t755*t758
1793 t762 = t322*t254
1794 t763 = t742 + t364 + t743 + t480 + t745 - 0.40e1_dp*t746 + t485 - 0.40e1_dp &
1795 *t748 - 0.1200e2_dp*t753 - 0.80e1_dp*t759 - 0.40e1_dp*t491 + t762
1796 t764 = t317*t254
1797 t766 = t729*t126
1798 t767 = t352*t327
1799 t769 = t732*t257
1800 t771 = t356*af
1801 t772 = t125*t133
1802 t773 = t772*t309
1803 t774 = t771*t773
1804 t777 = t143*af
1805 t778 = t125*t1
1806 t779 = t778*t337
1807 t780 = t777*t779
1808 t782 = t316*t305
1809 t783 = t305*t309
1810 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.40e1_dp*t784 + t705
1814
1815 e_ab(ip) = e_ab(ip) + ((t763 + t786)*t2 + t148 + t709 + t205 - t207 + t712 + t256 &
1816 + t259 + t715 + t321 + t323 - t325 + t326 + t329)*sc
1817 t789 = t309**2
1818 t793 = 0.2e1_dp*t131 + 0.2e1_dp*t338
1819 t796 = t312**2
1820 t799 = -t793
1821 t802 = 0.8549604655e0_dp*t332*t789 + 0.2564881399e1_dp*t129*t793 &
1822 + 0.8549604655e0_dp*t343*t796 + 0.2564881399e1_dp*t136*t799
1823 t803 = t802*t12
1824 t804 = t732*t309
1825 t806 = t356*t789
1826 t808 = t143*t793
1827 t811 = (-t803 - 0.80e1_dp*t804 - 0.1200e2_dp*t806 - 0.40e1_dp*t808)* &
1828 ap
1829 t812 = t811*t54
1830 t813 = 0.2e1_dp*t743
1831 t814 = t802*t56
1832 t815 = t814*t92
1833 t816 = t804*t92
1834 t817 = 0.80e1_dp*t816
1835 t818 = 0.2e1_dp*t762
1836 t819 = t806*t92
1837 t820 = 0.1200e2_dp*t819
1838 t821 = t808*t92
1839 t822 = 0.40e1_dp*t821
1840 t823 = 0.80e1_dp*t764
1841 t824 = t803*t126
1842 t825 = t732*t327
1843 t826 = 0.80e1_dp*t825
1844 t827 = 0.2e1_dp*t782
1845 t828 = t126*t789
1846 t829 = t356*t828
1847 t830 = 0.1200e2_dp*t829
1848 t831 = 0.80e1_dp*t784
1849 t832 = t126*t793
1850 t833 = t143*t832
1851 t834 = 0.40e1_dp*t833
1852 t835 = t812 + t813 + t480 + t815 - t817 + t818 - t820 - t822 - t823 &
1853 + 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 + t712 + 0.2e1_dp*t326 + 0.80e1_dp*t328 + t715)*sc
1857
1858 END IF
1859 IF (order >= 3 .OR. order == -3) THEN
1860 t842 = 0.3e1_dp*t480
1861 t851 = 0.3e1_dp*t592
1862 t857 = 0.3e1_dp*t705
1863 t859 = t17**(-0.2666666667e1_dp)
1864 t862 = 0.1e1_dp/t368/0.3141592654e1_dp
1865 t864 = 0.1e1_dp/t370/t130
1866 t865 = t862*t864
1867 t866 = t859*t24*t865
1868 t869 = t372*t164
1869 t870 = t366*t155*t869
1870 t873 = 0.1e1_dp/t370/t2
1871 t874 = t369*t873
1872 t875 = t367*t874
1873 t878 = t151*t385
1874 t881 = t379*t164
1875 t884 = t151*t399
1876 t887 = t16*t371
1877 t890 = t154**2
1878 t891 = 0.1e1_dp/t890
1879 t893 = t385*t164
1880 t896 = t164*t399
1881 t901 = 0.3365038134e0_dp*t859*t862*t864
1882 t903 = 0.1211413728e1_dp*t388*t873
1883 t905 = 0.1817120592e1_dp*t157*t371
1884 t906 = t17**(-0.2833333333e1_dp)
1885 t914 = -t901 + t903 - t905 - 0.2427089633e0_dp*bp*t906*t865 + 0.7943202439e0_dp &
1886 *t394*t874 - 0.9531842928e0_dp*t161*t887
1887 t937 = t17**(-0.2666666666e1_dp)
1888 t942 = t906*t862*t864
1889 t945 = t443*t873
1890 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 *t896 - t189*t914
1898 t972 = t39**(-0.40e1_dp)
1899 t979 = t874*t179
1900 t982 = 0.1e1_dp/t427
1901 t984 = t17**(-0.2500000000e1_dp)
1902 t986 = t865*t179
1903 t993 = 0.1e1_dp/t427/t172
1904 t996 = t865*t434
1905 t1010 = t887*t179
1906 t1013 = t874*t434
1907 t1019 = t427**2
1908 t1020 = 0.1e1_dp/t1019
1909 t1025 = t176**2
1910 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 *t45*t1020*t984*t1027
1922 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 *bp*t1020*t984*t1027
1939 t1085 = t15*t1084
1940 t1086 = t147*t479
1941 t1088 = t5**(-0.1666666667e1_dp)
1942 t1089 = t333*t133
1943 t1094 = t1*t371
1944 t1096 = 0.6e1_dp*t337 - 0.6e1_dp*t1094
1945 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 t1096
1950 t1109 = t1108*t12
1951 t1110 = t350*t142
1952 t1111 = t1110*t133
1953 t1113 = t140*t355
1954 t1114 = t1113*t333
1955 t1116 = t352*t340
1956 t1118 = t4**0.10e1_dp
1957 t1119 = t11*t1118
1958 t1120 = t1119*t1089
1959 t1125 = t143*t1096
1960 t1137 = t351*t305
1961 t1142 = t352*t601
1962 t1145 = t859*t97*t865
1963 t1148 = t372*t269
1964 t1149 = t366*t264*t1148
1965 t1151 = t607*t874
1966 t1154 = t151*t619
1967 t1157 = t379*t269
1968 t1160 = t151*t627
1969 t1165 = t263**2
1970 t1166 = 0.1e1_dp/t1165
1971 t1168 = t619*t269
1972 t1171 = t269*t627
1973 t1181 = -t901 + t903 - t905 - 0.2427089633e0_dp*bf*t906*t865 + 0.7943202439e0_dp &
1974 *t622*t874 - 0.9531842928e0_dp*t266*t887
1975 t1205 = t874*t283
1976 t1208 = 0.1e1_dp/t654
1977 t1211 = t865*t283
1978 t1218 = 0.1e1_dp/t654/t276
1979 t1221 = t865*t661
1980 t1235 = t887*t283
1981 t1238 = t874*t661
1982 t1244 = t654**2
1983 t1245 = 0.1e1_dp/t1244
1984 t1250 = t280**2
1985 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 *t674*t1171 - t291*t1181
1993 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 *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 *bf*t1245*t984*t1252 - t109*t1341*t122
2021 t1346 = t13*af*t1344
2022 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 *t750*t90*t91*t340 - 0.40e1_dp*t1125*t92 + 0.3600e2_dp*t1358
2030 t1361 = t362*t202
2031 t1363 = t353*t254
2032 t1365 = t144*t591
2033 t1368 = t859*t61*t865
2034 t1371 = t372*t217
2035 t1372 = t366*t212*t1371
2036 t1374 = t493*t874
2037 t1377 = t151*t505
2038 t1380 = t379*t217
2039 t1383 = t151*t513
2040 t1388 = t211**2
2041 t1389 = 0.1e1_dp/t1388
2042 t1391 = t505*t217
2043 t1394 = t217*t513
2044 t1404 = -t901 + t903 - t905 - 0.2427089633e0_dp*ba*t906*t865 + 0.7943202439e0_dp &
2045 *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 *t560*t1394 - t239*t1404
2053 t1469 = t74**(-0.40e1_dp)
2054 t1477 = t874*t231
2055 t1480 = 0.1e1_dp/t540
2056 t1483 = t865*t231
2057 t1490 = 0.1e1_dp/t540/t224
2058 t1493 = t865*t547
2059 t1507 = t887*t231
2060 t1510 = t874*t547
2061 t1516 = t540**2
2062 t1517 = 0.1e1_dp/t1516
2063 t1522 = t228**2
2064 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 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 t859*t1493 + 0.5132002385e0_dp*ba*t1517*t984*t1524
2092 t1570 = t57*aa*t1567*t91
2093 t1576 = t481*t254
2094 t1578 = t357*t254
2095 t1580 = t204*t591
2096 t1582 = t359*t254
2097 t1588 = t143*t305*t340
2098 t1590 = t141*t704
2099 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 *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 t1599)*t2)*sc
2111
2112 t1602 = 0.2400e2_dp*t753
2113 t1603 = 0.2e1_dp*t766
2114 t1604 = 0.80e1_dp*t746
2115 t1605 = 0.2e1_dp*t745
2116 t1606 = 0.80e1_dp*t748
2117 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 + t842 - t1609
2121 t1611 = 0.80e1_dp*t767
2122 t1613 = t322*t591
2123 t1615 = 0.80e1_dp*t732*t601
2124 t1620 = 0.80e1_dp*t733*t254
2125 t1622 = 0.80e1_dp*t352*t783
2126 t1626 = t732*t340
2127 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 *t136*t1639
2133 t1654 = t1653*t12
2134 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 t1142 - 0.40e1_dp*t1626*t92 + t1346 + t1654*t126 + 0.40e1_dp*t1657
2138 t1664 = 0.160e2_dp*t777*t304*t1*t337
2139 t1666 = 0.80e1_dp*t730*t254
2140 t1668 = t143*t1639
2141 t1671 = t728*t142
2142 t1672 = t1671*t133
2143 t1675 = t1110*t309
2144 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 *t1639
2150 t1696 = t316*t704
2151 t1697 = t315*t355
2152 t1698 = t1697*t333
2153 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.2e1_dp*t1580 - 0.40e1_dp*t1582 + 0.80e1_dp*t1671*t257
2159 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 t1626 - 0.1200e2_dp*t356*t340*t309 - 0.40e1_dp*t1668
2164 t1754 = 0.2e1_dp*t744*t254
2165 t1764 = 0.2e1_dp*t741*t202
2166 t1765 = t320*t479
2167 t1768 = 0.2400e2_dp*t750*t253*t751
2168 t1775 = 0.2e1_dp*t729*t305
2169 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 *t125*t340*t309 + t1775 - 0.40e1_dp*t1776
2174 t1782 = 0.2e1_dp*t742
2175 t1783 = 0.160e2_dp*t780
2176 t1785 = 0.2400e2_dp*t774
2177 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 - t484 + t1785 + t1786 + 0.4e1_dp*t485 + 0.160e2_dp*t602
2181
2182 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 *t743 - 0.160e2_dp*t764 - t492 + t365
2187 t1797 = t143*t305*t793
2188 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 *t725*t371
2195 t1829 = t352*t793
2196 t1832 = t814*t254
2197 t1836 = t732*t783
2198 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 *t1701*t772*t789 + 0.80e1_dp*t1836 + t1346 + t1838
2202 t1850 = t806*t254
2203 t1853 = t356*t305*t789
2204 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 *t133*t90*t91*t789 + 0.80e1_dp*t1671*t327
2210 t1870 = t804*t254
2211 t1872 = t143*t337
2212 t1884 = t1826*t12
2213 t1886 = t1671*t309
2214 t1891 = t808*t254
2215 t1893 = t803*t305
2216 t1894 = t1113*t789
2217 t1897 = t1858*t133
2218 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 *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 *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 *t1776
2234 t1949 = t842 - t1609 + t815 + t1611 + t851 + t1782 + t1783 + t597 + &
2235 t857 + (t1839 + t1868 + t1904 + t1945)*t2 - t817 + t1785 + t1786 &
2236 - t820 + t486 + t603 - t822
2237
2238 e_abb(ip) = e_abb(ip) + (t1795 + t1949)*sc
2239
2240 t1975 = t1697*t789
2241 t1979 = t789*t309
2242 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 *t136*t1986
2247 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 *t1979
2254 t2011 = t1858*t309
2255 t2014 = t732*t793
2256 t2019 = t143*t1986
2257 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.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 t2048)*t2)*sc
2271
2272 END IF
2273 END IF
2274 END DO
2275
2276!$OMP END DO
2277
2278 END SUBROUTINE vwn5_lsd_calc
2279
2280END MODULE xc_vwn
2281
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:28
bibliography::vosko1980
integer, save, public vosko1980
Definition bibliography.F:43
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition xc_derivative_desc.F:26
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:66
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_functionals_utilities
Utility routines for the functional calculations.
Definition xc_functionals_utilities.F:15
xc_functionals_utilities::set_util
subroutine, public set_util(cutoff)
...
Definition xc_functionals_utilities.F:54
xc_functionals_utilities::calc_srs_pw
subroutine, public calc_srs_pw(rho, x, n)
...
Definition xc_functionals_utilities.F:144
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::do_vwn5
integer, parameter, public do_vwn5
Definition xc_input_constants.F:30
xc_input_constants::do_vwn3
integer, parameter, public do_vwn3
Definition xc_input_constants.F:30
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition xc_rho_set_types.F:281
xc_vwn
Calculate the LDA functional according to Vosk, Wilk and Nusair Literature: S. H. Vosko,...
Definition xc_vwn.F:21
xc_vwn::vwn_lda_info
subroutine, public vwn_lda_info(reference, shortform, needs, max_deriv)
...
Definition xc_vwn.F:115
xc_vwn::vwn_lda_eval
subroutine, public vwn_lda_eval(rho_set, deriv_set, order, vwn_params)
...
Definition xc_vwn.F:167
xc_vwn::vwn_lsd_eval
subroutine, public vwn_lsd_eval(rho_set, deriv_set, order, vwn_params)
...
Definition xc_vwn.F:525
xc_vwn::vwn_lsd_info
subroutine, public vwn_lsd_info(reference, shortform, needs, max_deriv)
...
Definition xc_vwn.F:141
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:49
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78