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