(git:34ef472)
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 ! **************************************************************************************************
21 MODULE xc_vwn
22  USE bibliography, ONLY: vosko1980,&
23  cite_reference
24  USE input_section_types, ONLY: section_vals_type,&
25  section_vals_val_get
26  USE kinds, ONLY: dp
27  USE xc_derivative_desc, ONLY: deriv_rho,&
28  deriv_rhoa,&
29  deriv_rhob
30  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
31  xc_dset_get_derivative
32  USE xc_derivative_types, ONLY: xc_derivative_get,&
33  xc_derivative_type
34  USE xc_functionals_utilities, ONLY: calc_srs_pw,&
35  set_util
36  USE xc_input_constants, ONLY: do_vwn3,&
37  do_vwn5
38  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
39  USE xc_rho_set_types, ONLY: xc_rho_set_get,&
40  xc_rho_set_type
41 #include "../base/base_uses.f90"
42 
43  IMPLICIT NONE
44 
45  PRIVATE
46 
47 ! *** Global parameters ***
48 
49  REAL(KIND=dp), PARAMETER :: a = 0.0310907_dp, &
50  af = 0.01554535_dp, &
51  aa = -0.0168868639404_dp !-1/(6pi^2)
52 
53  PUBLIC :: vwn_lda_info, vwn_lda_eval, vwn_lsd_eval, vwn_lsd_info
54 
55  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_vwn'
56 
57  REAL(KIND=dp) :: eps_rho, b, c, x0, bf, cf, x0f, ba, ca, x0a
58 
59 CONTAINS
60 
61 ! **************************************************************************************************
62 !> \brief ...
63 !> \param cutoff ...
64 !> \param vwn_params ...
65 ! **************************************************************************************************
66  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 
2280 END 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