(git:e7e05ae)
xc_pbe.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 calculates the pbe correlation functional
10 !> \note
11 !> This was generated with the help of a maple worksheet that you can
12 !> find under doc/pbe.mw .
13 !> I did not add 3. derivatives in the polazied (lsd) case because it
14 !> would have added 2500 lines of code. If you need them the worksheet
15 !> is already prepared for them, and by uncommenting a couple of lines
16 !> you should be able to generate them.
17 !> \par History
18 !> 09.2004 created [fawzi]
19 !> \author fawzi
20 ! **************************************************************************************************
21 MODULE xc_pbe
22  USE bibliography, ONLY: perdew1996,&
23  perdew2008,&
24  zhang1998,&
25  cite_reference
26  USE input_section_types, ONLY: section_vals_type,&
27  section_vals_val_get
28  USE kinds, ONLY: dp
29  USE mathconstants, ONLY: pi
30  USE xc_derivative_desc, ONLY: deriv_norm_drho,&
31  deriv_norm_drhoa,&
32  deriv_norm_drhob,&
33  deriv_rho,&
34  deriv_rhoa,&
35  deriv_rhob
36  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
37  xc_dset_get_derivative
38  USE xc_derivative_types, ONLY: xc_derivative_get,&
39  xc_derivative_type
40  USE xc_input_constants, ONLY: xc_pbe_orig,&
41  xc_pbe_rev,&
42  xc_pbe_sol
43  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
44  USE xc_rho_set_types, ONLY: xc_rho_set_get,&
45  xc_rho_set_type
46 #include "../base/base_uses.f90"
47 
48  IMPLICIT NONE
49  PRIVATE
50 
51  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
52  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pbe'
53  REAL(kind=dp), PARAMETER, PRIVATE :: a = 0.04918_dp, b = 0.132_dp, &
54  c = 0.2533_dp, d = 0.349_dp
55 
56  PUBLIC :: pbe_lda_info, pbe_lsd_info, pbe_lda_eval, pbe_lsd_eval
57 CONTAINS
58 
59 ! **************************************************************************************************
60 !> \brief return various information on the functional
61 !> \param pbe_params section selecting the various parameters for the functional
62 !> \param reference string with the reference of the actual functional
63 !> \param shortform string with the shortform of the functional name
64 !> \param needs the components needed by this functional are set to
65 !> true (does not set the unneeded components to false)
66 !> \param max_deriv ...
67 !> \author fawzi
68 ! **************************************************************************************************
69  SUBROUTINE pbe_lda_info(pbe_params, reference, shortform, needs, max_deriv)
70  TYPE(section_vals_type), POINTER :: pbe_params
71  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
72  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
73  INTEGER, INTENT(out), OPTIONAL :: max_deriv
74 
75  INTEGER :: param
76  REAL(kind=dp) :: sc, sx
77 
78  CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
79  CALL section_vals_val_get(pbe_params, "scale_x", r_val=sx)
80  CALL section_vals_val_get(pbe_params, "scale_c", r_val=sc)
81 
82  SELECT CASE (param)
83  CASE (xc_pbe_orig)
84  CALL cite_reference(perdew1996)
85  IF (sx == 1._dp .AND. sc == 1._dp) THEN
86  IF (PRESENT(reference)) THEN
87  reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "// &
88  "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996) {spin unpolarized}"
89  END IF
90  IF (PRESENT(shortform)) THEN
91  shortform = "PBE Perdew-Burke-Ernzerhof xc functional (unpolarized)"
92  END IF
93  ELSE
94  IF (PRESENT(reference)) THEN
95  WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
96  "J.P.Perdew, K.Burke, M.Ernzerhof, ", &
97  "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)", &
98  sx, sc, "{spin unpolarized}"
99  END IF
100  IF (PRESENT(shortform)) THEN
101  WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") &
102  "PBE, Perdew-Burke-Ernzerhof xc functional (unpolarized)", sx, sc
103  END IF
104  END IF
105  CASE (xc_pbe_rev)
106  CALL cite_reference(perdew1996)
107  CALL cite_reference(zhang1998)
108  IF (sx == 1._dp .AND. sc == 1._dp) THEN
109  IF (PRESENT(reference)) THEN
110  reference = "revPBE, Yingkay Zhang and Weitao Yang,"// &
111  " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998) {spin unpolarized}"
112  END IF
113  IF (PRESENT(shortform)) THEN
114  shortform = "revPBE, revised PBE xc functional (unpolarized)"
115  END IF
116  ELSE
117  IF (PRESENT(reference)) THEN
118  WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
119  "revPBE, Yingkay Zhang and Weitao Yang,", &
120  " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)", &
121  sx, sc, "{spin unpolarized}"
122  END IF
123  IF (PRESENT(shortform)) THEN
124  WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") &
125  "revPBE, revised PBE xc functional (unpolarized)", sx, sc
126  END IF
127  END IF
128  CASE (xc_pbe_sol)
129  CALL cite_reference(perdew1996)
130  CALL cite_reference(perdew2008)
131  IF (sx == 1._dp .AND. sc == 1._dp) THEN
132  IF (PRESENT(reference)) THEN
133  reference = "PBEsol, J.P. Perdew et al., "// &
134  "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) "// &
135  "{spin unpolarized}"
136  END IF
137  IF (PRESENT(shortform)) THEN
138  shortform = "PBEsol, PBE for solids and surfaces xc functional (unpolarized)"
139  END IF
140  ELSE
141  IF (PRESENT(reference)) THEN
142  WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
143  "PBEsol, J.P. Perdew et al., ", &
144  "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) ", &
145  sx, sc, "{spin unpolarized}"
146  END IF
147  IF (PRESENT(shortform)) THEN
148  WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") &
149  "PBEsol, PBE for solids and surfaces xc functional (unpolarized)", sx, sc
150  END IF
151  END IF
152  CASE default
153  cpabort("Unsupported parametrization")
154  END SELECT
155  IF (PRESENT(needs)) THEN
156  needs%rho = .true.
157  needs%norm_drho = .true.
158  END IF
159  IF (PRESENT(max_deriv)) max_deriv = 3
160 
161  END SUBROUTINE pbe_lda_info
162 
163 ! **************************************************************************************************
164 !> \brief return various information on the functional
165 !> \param pbe_params section selecting the various parameters for the functional
166 !> \param reference string with the reference of the actual functional
167 !> \param shortform string with the shortform of the functional name
168 !> \param needs the components needed by this functional are set to
169 !> true (does not set the unneeded components to false)
170 !> \param max_deriv ...
171 !> \author fawzi
172 ! **************************************************************************************************
173  SUBROUTINE pbe_lsd_info(pbe_params, reference, shortform, needs, max_deriv)
174  TYPE(section_vals_type), POINTER :: pbe_params
175  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
176  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
177  INTEGER, INTENT(out), OPTIONAL :: max_deriv
178 
179  INTEGER :: param
180  REAL(kind=dp) :: sc, sx
181 
182  CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
183  CALL section_vals_val_get(pbe_params, "scale_x", r_val=sx)
184  CALL section_vals_val_get(pbe_params, "scale_c", r_val=sc)
185 
186  SELECT CASE (param)
187  CASE (xc_pbe_orig)
188  CALL cite_reference(perdew1996)
189  IF (sx == 1._dp .AND. sc == 1._dp) THEN
190  IF (PRESENT(reference)) THEN
191  reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "// &
192  "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996) {spin polarized}"
193  END IF
194  IF (PRESENT(shortform)) THEN
195  shortform = "PBE xc functional (spin polarized)"
196  END IF
197  ELSE
198  IF (PRESENT(reference)) THEN
199  WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
200  "J.P.Perdew, K.Burke, M.Ernzerhof, ", &
201  "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)", &
202  sx, sc, "{spin polarized}"
203  END IF
204  IF (PRESENT(shortform)) THEN
205  WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") &
206  "PBE xc functional (spin polarized)", sx, sc
207  END IF
208  END IF
209  CASE (xc_pbe_rev)
210  CALL cite_reference(perdew1996)
211  CALL cite_reference(zhang1998)
212  IF (sx == 1._dp .AND. sc == 1._dp) THEN
213  IF (PRESENT(reference)) THEN
214  reference = "revPBE, Yingkay Zhang and Weitao Yang,"// &
215  " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998) {spin polarized}"
216  END IF
217  IF (PRESENT(shortform)) THEN
218  shortform = "revPBE, revised PBE xc functional (spin polarized)"
219  END IF
220  ELSE
221  IF (PRESENT(reference)) THEN
222  WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
223  "revPBE, Yingkay Zhang and Weitao Yang,", &
224  " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)", &
225  sx, sc, "{spin polarized}"
226  END IF
227  IF (PRESENT(shortform)) THEN
228  WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") &
229  "revPBE, revised PBE xc functional (spin polarized)", sx, sc
230  END IF
231  END IF
232  CASE (xc_pbe_sol)
233  CALL cite_reference(perdew1996)
234  CALL cite_reference(perdew2008)
235  IF (sx == 1._dp .AND. sc == 1._dp) THEN
236  IF (PRESENT(reference)) THEN
237  reference = "PBEsol, J.P. Perdew et al., "// &
238  "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) "// &
239  "{spin polarized}"
240  END IF
241  IF (PRESENT(shortform)) THEN
242  shortform = "PBEsol, PBE for solids and surfaces xc functional (spin polarized)"
243  END IF
244  ELSE
245  IF (PRESENT(reference)) THEN
246  WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
247  "PBEsol, J.P. Perdew et al., ", &
248  "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) ", &
249  sx, sc, "{spin unpolarized}"
250  END IF
251  IF (PRESENT(shortform)) THEN
252  WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") &
253  "PBEsol, PBE for solids and surfaces xc functional (spin polarized)", sx, sc
254  END IF
255  END IF
256  CASE default
257  cpabort("Unsupported parametrization")
258  END SELECT
259  IF (PRESENT(needs)) THEN
260  needs%rho_spin = .true.
261  needs%norm_drho_spin = .true.
262  needs%norm_drho = .true.
263  END IF
264  IF (PRESENT(max_deriv)) max_deriv = 2
265 
266  END SUBROUTINE pbe_lsd_info
267 
268 ! **************************************************************************************************
269 !> \brief evaluates the pbe correlation functional for lda
270 !> \param rho_set the density where you want to evaluate the functional
271 !> \param deriv_set place where to store the functional derivatives (they are
272 !> added to the derivatives)
273 !> \param grad_deriv degree of the derivative that should be evaluated,
274 !> if positive all the derivatives up to the given degree are evaluated,
275 !> if negative only the given degree is calculated
276 !> \param pbe_params ...
277 !> \author fawzi
278 ! **************************************************************************************************
279  SUBROUTINE pbe_lda_eval(rho_set, deriv_set, grad_deriv, pbe_params)
280  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
281  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
282  INTEGER, INTENT(in) :: grad_deriv
283  TYPE(section_vals_type), POINTER :: pbe_params
284 
285  CHARACTER(len=*), PARAMETER :: routinen = 'pbe_lda_eval'
286 
287  INTEGER :: handle, npoints, param
288  INTEGER, DIMENSION(2, 3) :: bo
289  REAL(kind=dp) :: epsilon_rho, scale_ec, scale_ex
290  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
291  e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
292  e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
293  TYPE(xc_derivative_type), POINTER :: deriv
294 
295  CALL timeset(routinen, handle)
296 
297  CALL xc_rho_set_get(rho_set, rho=rho, &
298  norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
299  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
300 
301  dummy => rho
302 
303  e_0 => dummy
304  e_rho => dummy
305  e_ndrho => dummy
306  e_rho_rho => dummy
307  e_ndrho_rho => dummy
308  e_ndrho_ndrho => dummy
309  e_rho_rho_rho => dummy
310  e_ndrho_rho_rho => dummy
311  e_ndrho_ndrho_rho => dummy
312  e_ndrho_ndrho_ndrho => dummy
313 
314  IF (grad_deriv >= 0) THEN
315  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
316  allocate_deriv=.true.)
317  CALL xc_derivative_get(deriv, deriv_data=e_0)
318  END IF
319  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
320  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
321  allocate_deriv=.true.)
322  CALL xc_derivative_get(deriv, deriv_data=e_rho)
323  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
324  allocate_deriv=.true.)
325  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
326  END IF
327  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
328  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
329  allocate_deriv=.true.)
330  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
331  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
332  allocate_deriv=.true.)
333  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
334  deriv => xc_dset_get_derivative(deriv_set, &
335  [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
336  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
337  END IF
338  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
339  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
340  allocate_deriv=.true.)
341  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
342  deriv => xc_dset_get_derivative(deriv_set, &
343  [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.true.)
344  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
345  deriv => xc_dset_get_derivative(deriv_set, &
346  [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.true.)
347  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
348  deriv => xc_dset_get_derivative(deriv_set, &
349  [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
350  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
351  END IF
352  IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
353  cpabort("derivatives bigger than 3 not implemented")
354  END IF
355 
356  CALL section_vals_val_get(pbe_params, "scale_c", r_val=scale_ec)
357  CALL section_vals_val_get(pbe_params, "scale_x", r_val=scale_ex)
358  CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
359 
360 !$OMP PARALLEL DEFAULT(NONE), &
361 !$OMP SHARED(rho,norm_drho,e_0,e_rho,e_ndrho,e_rho_rho,e_ndrho_rho),&
362 !$OMP SHARED(e_ndrho_ndrho,e_rho_rho_rho,e_ndrho_rho_rho,e_ndrho_ndrho_rho),&
363 !$OMP SHARED(e_ndrho_ndrho_ndrho,grad_deriv,npoints,epsilon_rho),&
364 !$OMP SHARED(pbe_params),&
365 !$OMP SHARED(param,scale_ec,scale_ex)
366  CALL pbe_lda_calc(rho=rho, norm_drho=norm_drho, &
367  e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
368  e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
369  e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
370  e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, &
371  grad_deriv=grad_deriv, &
372  npoints=npoints, epsilon_rho=epsilon_rho, &
373  param=param, scale_ec=scale_ec, scale_ex=scale_ex)
374 !$OMP END PARALLEL
375 
376  CALL timestop(handle)
377  END SUBROUTINE pbe_lda_eval
378 
379 ! **************************************************************************************************
380 !> \brief evaluates the pbe correlation functional for lda
381 !> \param rho the density where you want to evaluate the functional
382 !> \param norm_drho ...
383 !> \param e_0 ...
384 !> \param e_rho ...
385 !> \param e_ndrho ...
386 !> \param e_rho_rho ...
387 !> \param e_ndrho_rho ...
388 !> \param e_ndrho_ndrho ...
389 !> \param e_rho_rho_rho ...
390 !> \param e_ndrho_rho_rho ...
391 !> \param e_ndrho_ndrho_rho ...
392 !> \param e_ndrho_ndrho_ndrho ...
393 !> \param grad_deriv degree of the derivative that should be evaluated,
394 !> if positive all the derivatives up to the given degree are evaluated,
395 !> if negative only the given degree is calculated
396 !> \param npoints ...
397 !> \param epsilon_rho ...
398 !> \param ! ...
399 !> \param pbe_params parameters for the pbe functional
400 !> \author fawzi
401 ! **************************************************************************************************
402  SUBROUTINE pbe_lda_calc(rho, norm_drho, &
403  e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
404  e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
405  e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, &
406  ! pbe_params)
407  param, scale_ec, scale_ex)
408  INTEGER, INTENT(in) :: npoints, grad_deriv
409  REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_ndrho, &
410  e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, &
411  e_ndrho, e_rho, e_0
412  REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho
413  REAL(kind=dp), INTENT(in) :: epsilon_rho
414  INTEGER, INTENT(in) :: param
415  REAL(kind=dp), INTENT(in) :: scale_ec, scale_ex
416 
417  INTEGER :: ii
418  REAL(kind=dp) :: a, a1rho, a1rhorho, a2rho, a_1, alpha_1_1, arho, arho1rho, arhorho, &
419  arhorhorho, beta, beta_1_1, beta_2_1, beta_3_1, beta_4_1, e_c_u_0, e_c_u_01rho, &
420  e_c_u_01rhorho, e_c_u_02rho, e_c_u_0rho, e_c_u_0rho1rho, e_c_u_0rhorho, e_c_u_0rhorhorho, &
421  epsilon_cgga, epsilon_cggarho, epsilon_cggarhorho, ex_ldarhorhorho, ex_unif, ex_unif1rho, &
422  ex_unif1rhorho, ex_unif2rho, ex_unifrho, ex_unifrho1rho, ex_unifrhorho, fx, fx1rho, &
423  fx1rhorho, fx2rho, fxnorm_drho, fxnorm_drho1rho, fxnorm_drhonorm_drho, fxnorm_drhorho, &
424  fxrho, fxrho1rho, fxrhorho, gamma_var, hnorm_drho, hnorm_drhonorm_drho
425  REAL(kind=dp) :: hnorm_drhorho, k_f, k_f2rho, k_frho, k_frhorho, k_frhorhorho, k_s, k_s1rho, &
426  k_s1rhorho, k_s2rho, k_srho, k_srho1rho, k_srhorho, kappa, kf, kf2rho, kfrho, kfrhorho, &
427  kfrhorhorho, mu, my_norm_drho, my_rho, p_1, p_2, rs, rs2rho, rsrho, rsrhorho, &
428  rsrhorhorho, s, s1rho, s1rhorho, s2norm_drho, s2rho, snorm_drho, snorm_drho1rho, &
429  snorm_drhorho, srho, srho1rho, srhorho, t, t1, t1001, t1004, t1005, t1006, t1008, t101, &
430  t1011, t1012, t1013, t1014, t1016, t1017, t1018, t1019, t1022, t1024, t1026, t1028, t103, &
431  t1030, t1031, t1035, t1042, t1046, t1048, t105, t1050, t1052, t1054, t1055
432  REAL(kind=dp) :: t1060, t1061, t1062, t1067, t108, t1093, t1097, t11, t1103, t1104, t1106, &
433  t1109, t111, t1115, t1118, t1121, t1122, t1126, t1129, t113, t114, t1148, t115, t1152, &
434  t1157, t1167, t1187, t119, t1196, t1203, t1218, t1226, t123, t124, t1249, t125, t126, &
435  t1262, t1263, t127, t1291, t1292, t1295, t13, t131, t1329, t1342, t135, t138, t1380, &
436  t1385, t1386, t1387, t1389, t139, t1393, t14, t140, t142, t146, t1468, t148, t1482, t149, &
437  t1492, t1493, t150, t1504, t1505, t151, t1511, t1515, t1521, t1525, t1528, t153, t1532, &
438  t1545, t155, t1565, t157, t1573, t158, t1584, t159, t1608, t1612
439  REAL(kind=dp) :: t162, t1628, t1629, t163, t1633, t164, t1646, t1652, t1660, t167, t1672, &
440  t1676, t168, t17, t170, t171, t1715, t1722, t175, t1758, t1759, t176, t1761, t177, t178, &
441  t179, t1797, t182, t183, t1838, t1851, t186, t187, t1878, t1889, t189, t19, t190, t1907, &
442  t191, t1922, t1927, t1933, t1937, t195, t1952, t196, t1964, t1965, t1968, t1969, t197, &
443  t1972, t198, t1990, t1rho, t1rhorho, t2, t20, t200, t202, t2020, t2024, t2028, t2031, &
444  t204, t2041, t208, t21, t214, t217, t218, t22, t226, t229, t230, t238, t239, t240, t241, &
445  t246, t252, t253, t255, t259, t26, t260, t261, t262, t265, t266
446  REAL(kind=dp) :: t267, t27, t271, t272, t273, t277, t278, t280, t281, t286, t290, t291, &
447  t293, t294, t295, t296, t297, t299, t2norm_drho, t2rho, t3, t305, t309, t310, t315, t317, &
448  t318, t321, t323, t324, t327, t329, t330, t331, t336, t338, t339, t340, t341, t348, t349, &
449  t354, t357, t359, t360, t361, t362, t369, t370, t371, t374, t377, t378, t381, t382, t384, &
450  t385, t387, t388, t390, t392, t393, t397, t4, t400, t401, t404, t408, t409, t410, t414, &
451  t417, t418, t423, t426, t427, t432, t435, t436, t438, t439, t440, t448, t449, t451, t456, &
452  t457, t458, t459, t461, t462, t463, t465, t466, t469, t470
453  REAL(kind=dp) :: t471, t472, t476, t487, t491, t496, t498, t5, t500, t503, t506, t507, t510, &
454  t513, t517, t521, t525, t529, t530, t535, t541, t548, t553, t556, t557, t559, t562, t566, &
455  t581, t586, t590, t591, t594, t598, t6, t605, t609, t611, t612, t614, t627, t645, t65, &
456  t654, t656, t661, t664, t669, t670, t671, t673, t675, t685, t69, t693, t7, t70, t701, &
457  t702, t71, t714, t717, t72, t720, t73, t74, t740, t743, t748, t75, t758, t77, t776, t78, &
458  t8, t80, t801, t809, t81, t812, t821, t825, t83, t831, t84, t85, t86, t868, t87, t877, &
459  t879, t88, t880, t885, t90, t91, t94, t940, t942, t943, t945
460  REAL(kind=dp) :: t946, t948, t95, t950, t951, t954, t967, t976, t98, t980, t982, t984, t989, &
461  t99, t990, t994, t995, t998, t999, tnorm_drho, tnorm_drho1rho, tnorm_drhorho, &
462  tnorm_drhorhorho, trho, trho1rho, trhorho, trhorhorho
463 
464 !TYPE(section_vals_type), POINTER :: pbe_params
465 !, param
466 ! scale_ec, scale_ex, snorm_drho, snorm_drho1rho, snorm_drhorho, srho, &
467 ! Parametrization of PBE
468 
469  SELECT CASE (param)
470  ! Original PBE
471  CASE (xc_pbe_orig)
472  kappa = 0.804e0_dp
473  beta = 0.66725e-1_dp
474  mu = beta*pi**2/0.3e1_dp
475  ! Revised PBE (revPBE)
476  CASE (xc_pbe_rev)
477  kappa = 0.1245e1_dp
478  beta = 0.66725e-1_dp
479  mu = beta*pi**2/0.3e1_dp
480  ! PBE for solids (PBEsol)
481  CASE (xc_pbe_sol)
482  kappa = 0.804e0_dp
483  beta = 0.46e-1_dp
484  mu = 0.1e2_dp/0.81e2_dp
485  CASE default
486  cpabort("Unsupported parametrization")
487  END SELECT
488 
489  gamma_var = (0.1e1_dp - log(0.2e1_dp))/pi**2
490  p_1 = 0.10e1_dp
491  a_1 = 0.31091e-1_dp
492  alpha_1_1 = 0.21370e0_dp
493  beta_1_1 = 0.75957e1_dp
494  beta_2_1 = 0.35876e1_dp
495  beta_3_1 = 0.16382e1_dp
496  beta_4_1 = 0.49294e0_dp
497  p_2 = 0.10e1_dp
498  t1 = 3**(0.1e1_dp/0.3e1_dp)
499  t2 = 4**(0.1e1_dp/0.3e1_dp)
500  t3 = t2**2
501  t4 = t1*t3
502  t5 = 0.1e1_dp/pi
503  t69 = pi**2
504  t627 = 0.1e1_dp/t69/pi
505 
506  SELECT CASE (grad_deriv)
507  CASE default
508 !$OMP DO
509  DO ii = 1, npoints
510  my_rho = rho(ii)
511  IF (my_rho > epsilon_rho) THEN
512  my_norm_drho = norm_drho(ii)
513 
514  t6 = 0.1e1_dp/my_rho
515  t7 = t5*t6
516  t8 = t7**(0.1e1_dp/0.3e1_dp)
517  rs = t4*t8/0.4e1_dp
518  t11 = 0.1e1_dp + alpha_1_1*rs
519  t13 = 0.1e1_dp/a_1
520  t14 = sqrt(rs)
521  t17 = t14*rs
522  t19 = p_1 + 0.1e1_dp
523  t20 = rs**t19
524  t21 = beta_4_1*t20
525  t22 = beta_1_1*t14 + beta_2_1*rs + beta_3_1*t17 + t21
526  t26 = 0.1e1_dp + t13/t22/0.2e1_dp
527  t27 = log(t26)
528  e_c_u_0 = -0.2e1_dp*a_1*t11*t27
529  t65 = 2**(0.1e1_dp/0.3e1_dp)
530  t70 = t69*my_rho
531  t71 = t70**(0.1e1_dp/0.3e1_dp)
532  k_f = t1*t71
533  t72 = k_f*t5
534  t73 = sqrt(t72)
535  k_s = 0.2e1_dp*t73
536  t74 = 0.1e1_dp/k_s
537  t75 = my_norm_drho*t74
538  t = t75*t6/0.2e1_dp
539  t77 = 0.1e1_dp/gamma_var
540  t78 = beta*t77
541  t80 = exp(-e_c_u_0*t77)
542  t81 = -0.1e1_dp + t80
543  a = t78/t81
544  t83 = t**2
545  t84 = a*t83
546  t85 = 0.1e1_dp + t84
547  t86 = t83*t85
548  t87 = a**2
549  t88 = t83**2
550  t90 = 0.1e1_dp + t84 + t87*t88
551  t91 = 0.1e1_dp/t90
552  t94 = 0.1e1_dp + t78*t86*t91
553  t95 = log(t94)
554  epsilon_cgga = e_c_u_0 + gamma_var*t95
555  kf = k_f
556  ex_unif = -0.3e1_dp/0.4e1_dp*t5*kf
557  t98 = 0.1e1_dp/kf
558  t99 = my_norm_drho*t98
559  s = t99*t6/0.2e1_dp
560  t101 = s**2
561  t103 = 0.1e1_dp/kappa
562  t105 = 0.1e1_dp + mu*t101*t103
563  fx = 0.1e1_dp + kappa - kappa/t105
564  t108 = my_rho*ex_unif
565 
566  IF (grad_deriv >= 0) THEN
567  e_0(ii) = e_0(ii) + &
568  scale_ex*t108*fx + scale_ec*my_rho*epsilon_cgga
569  END IF
570 
571  t111 = t8**2
572  t113 = 0.1e1_dp/t111*t5
573  t114 = my_rho**2
574  t115 = 0.1e1_dp/t114
575  rsrho = -t4*t113*t115/0.12e2_dp
576  t119 = a_1*alpha_1_1
577  t123 = t22**2
578  t124 = 0.1e1_dp/t123
579  t125 = t11*t124
580  t126 = 0.1e1_dp/t14
581  t127 = beta_1_1*t126
582  t131 = beta_3_1*t14
583  t135 = 0.1e1_dp/rs
584  t138 = t127*rsrho/0.2e1_dp + beta_2_1*rsrho + 0.3e1_dp/ &
585  0.2e1_dp*t131*rsrho + t21*t19*rsrho*t135
586  t139 = 0.1e1_dp/t26
587  t140 = t138*t139
588  e_c_u_0rho = -0.2e1_dp*t119*rsrho*t27 + t125*t140
589  t142 = t71**2
590  k_frho = t1/t142*t69/0.3e1_dp
591  t146 = 0.1e1_dp/t73
592  k_srho = t146*k_frho*t5
593  t148 = k_s**2
594  t149 = 0.1e1_dp/t148
595  t150 = my_norm_drho*t149
596  t151 = t6*k_srho
597  t153 = t75*t115
598  trho = -t150*t151/0.2e1_dp - t153/0.2e1_dp
599  t155 = gamma_var**2
600  t157 = beta/t155
601  t158 = t81**2
602  t159 = 0.1e1_dp/t158
603  arho = t157*t159*e_c_u_0rho*t80
604  t162 = t78*t
605  t163 = t85*t91
606  t164 = t163*trho
607  t167 = arho*t83
608  t168 = a*t
609  t170 = 0.2e1_dp*t168*trho
610  t171 = t167 + t170
611  t175 = t78*t83
612  t176 = t90**2
613  t177 = 0.1e1_dp/t176
614  t178 = t85*t177
615  t179 = a*t88
616  t182 = t83*t
617  t183 = t87*t182
618  t186 = t167 + t170 + 0.2e1_dp*t179*arho + 0.4e1_dp*t183*trho
619  t187 = t178*t186
620  t189 = 0.2e1_dp*t162*t164 + t78*t83*t171*t91 - t175*t187
621  t190 = gamma_var*t189
622  t191 = 0.1e1_dp/t94
623  epsilon_cggarho = e_c_u_0rho + t190*t191
624  kfrho = k_frho
625  ex_unifrho = -0.3e1_dp/0.4e1_dp*t5*kfrho
626  t195 = kf**2
627  t196 = 0.1e1_dp/t195
628  t197 = my_norm_drho*t196
629  t198 = t6*kfrho
630  t200 = t99*t115
631  srho = -t197*t198/0.2e1_dp - t200/0.2e1_dp
632  t202 = t105**2
633  t204 = 0.1e1_dp/t202*mu
634  fxrho = 0.2e1_dp*t204*s*srho
635  t208 = my_rho*ex_unifrho
636 
637  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
638  e_rho(ii) = e_rho(ii) + &
639  scale_ex*(ex_unif*fx + t208*fx + t108*fxrho) + &
640  scale_ec*(epsilon_cgga + my_rho*epsilon_cggarho)
641  END IF
642 
643  tnorm_drho = t74*t6/0.2e1_dp
644  t214 = t163*tnorm_drho
645  t217 = t78*t182
646  t218 = a*tnorm_drho
647  t226 = 0.2e1_dp*t168*tnorm_drho + 0.4e1_dp*t183*tnorm_drho
648  t229 = 0.2e1_dp*t162*t214 + 0.2e1_dp*t217*t218*t91 - &
649  t175*t178*t226
650  t230 = gamma_var*t229
651  hnorm_drho = t230*t191
652  snorm_drho = t98*t6/0.2e1_dp
653  fxnorm_drho = 0.2e1_dp*t204*s*snorm_drho
654 
655  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
656  e_ndrho(ii) = e_ndrho(ii) + &
657  scale_ex*t108*fxnorm_drho + scale_ec*my_rho* &
658  hnorm_drho
659  END IF
660 
661  t238 = 0.1e1_dp/t69
662  t239 = 0.1e1_dp/t111/t7*t238
663  t240 = t114**2
664  t241 = 0.1e1_dp/t240
665  t246 = 0.1e1_dp/t114/my_rho
666  rsrhorho = -t4*t239*t241/0.18e2_dp + t4*t113* &
667  t246/0.6e1_dp
668  t252 = 0.2e1_dp*t119*rsrhorho*t27
669  t253 = alpha_1_1*rsrho
670  t255 = t124*t138*t139
671  t259 = 0.1e1_dp/t123/t22
672  t260 = t11*t259
673  t261 = t138**2
674  t262 = t261*t139
675  t265 = 0.1e1_dp/t17
676  t266 = beta_1_1*t265
677  t267 = rsrho**2
678  t271 = t127*rsrhorho/0.2e1_dp
679  t272 = beta_2_1*rsrhorho
680  t273 = beta_3_1*t126
681  t277 = 0.3e1_dp/0.2e1_dp*t131*rsrhorho
682  t278 = t19**2
683  t280 = rs**2
684  t281 = 0.1e1_dp/t280
685  t286 = t21*t19*rsrhorho*t135
686  t290 = -t266*t267/0.4e1_dp + t271 + t272 + 0.3e1_dp/0.4e1_dp &
687  *t273*t267 + t277 + t21*t278*t267*t281 + t286 - t21*t19 &
688  *t267*t281
689  t291 = t290*t139
690  t293 = t123**2
691  t294 = 0.1e1_dp/t293
692  t295 = t11*t294
693  t296 = t26**2
694  t297 = 0.1e1_dp/t296
695  t299 = t261*t297*t13
696  e_c_u_0rhorho = -t252 + 0.2e1_dp*t253*t255 - 0.2e1_dp*t260* &
697  t262 + t125*t291 + t295*t299/0.2e1_dp
698  e_c_u_01rho = e_c_u_0rho
699  t305 = t69**2
700  k_frhorho = -0.2e1_dp/0.9e1_dp*t1/t142/t70*t305
701  t309 = 0.1e1_dp/t73/t72
702  t310 = k_frho**2
703  t315 = t146*k_frhorho*t5
704  k_srhorho = -t309*t310*t238/0.2e1_dp + t315
705  k_s1rho = k_srho
706  t317 = 0.1e1_dp/t148/k_s
707  t318 = my_norm_drho*t317
708  t321 = t115*k_srho
709  t323 = t150*t321/0.2e1_dp
710  t324 = t6*k_srhorho
711  t327 = t115*k_s1rho
712  t329 = t150*t327/0.2e1_dp
713  t330 = t75*t246
714  trhorho = t318*t151*k_s1rho + t323 - t150*t324/0.2e1_dp + &
715  t329 + t330
716  t331 = t6*k_s1rho
717  t1rho = -t150*t331/0.2e1_dp - t153/0.2e1_dp
718  t336 = beta/t155/gamma_var
719  t338 = 0.1e1_dp/t158/t81
720  t339 = t336*t338
721  t340 = t80**2
722  t341 = e_c_u_0rho*t340
723  t348 = t336*t159
724  t349 = e_c_u_0rho*e_c_u_01rho
725  arhorho = 0.2e1_dp*t339*t341*e_c_u_01rho + t157*t159* &
726  e_c_u_0rhorho*t80 - t348*t349*t80
727  a1rho = t157*t159*e_c_u_01rho*t80
728  t354 = t78*t1rho
729  t357 = a1rho*t83
730  t359 = 0.2e1_dp*t168*t1rho
731  t360 = t357 + t359
732  t361 = t360*t91
733  t362 = t361*trho
734  t369 = t357 + t359 + 0.2e1_dp*t179*a1rho + 0.4e1_dp*t183*t1rho
735  t370 = trho*t369
736  t371 = t178*t370
737  t374 = t163*trhorho
738  t377 = t171*t91
739  t378 = t377*t1rho
740  t381 = arhorho*t83
741  t382 = arho*t
742  t384 = 0.2e1_dp*t382*t1rho
743  t385 = a1rho*t
744  t387 = 0.2e1_dp*t385*trho
745  t388 = a*t1rho
746  t390 = 0.2e1_dp*t388*trho
747  t392 = 0.2e1_dp*t168*trhorho
748  t393 = t381 + t384 + t387 + t390 + t392
749  t397 = t171*t177
750  t400 = t186*t1rho
751  t401 = t178*t400
752  t404 = t360*t177
753  t408 = 0.1e1_dp/t176/t90
754  t409 = t85*t408
755  t410 = t186*t369
756  t414 = a1rho*t88
757  t417 = a*t182
758  t418 = arho*t1rho
759  t423 = trho*a1rho
760  t426 = t87*t83
761  t427 = trho*t1rho
762  t432 = t381 + t384 + t387 + t390 + t392 + 0.2e1_dp*t414*arho + &
763  0.8e1_dp*t417*t418 + 0.2e1_dp*t179*arhorho + 0.8e1_dp* &
764  t417*t423 + 0.12e2_dp*t426*t427 + 0.4e1_dp*t183*trhorho
765  t435 = 0.2e1_dp*t354*t164 + 0.2e1_dp*t162*t362 - 0.2e1_dp &
766  *t162*t371 + 0.2e1_dp*t162*t374 + 0.2e1_dp*t162*t378 + &
767  t78*t83*t393*t91 - t175*t397*t369 - 0.2e1_dp*t162*t401 &
768  - t175*t404*t186 + 0.2e1_dp*t175*t409*t410 - t175*t178 &
769  *t432
770  t436 = gamma_var*t435
771  t438 = t94**2
772  t439 = 0.1e1_dp/t438
773  t440 = t163*t1rho
774  t448 = 0.2e1_dp*t162*t440 + t78*t83*t360*t91 - t175* &
775  t178*t369
776  t449 = t439*t448
777  t451 = gamma_var*t448
778  epsilon_cggarhorho = e_c_u_0rhorho + t436*t191 - t190*t449
779  kfrhorho = k_frhorho
780  ex_unifrhorho = -0.3e1_dp/0.4e1_dp*t5*kfrhorho
781  ex_unif1rho = ex_unifrho
782  t456 = 0.1e1_dp/t195/kf
783  t457 = my_norm_drho*t456
784  t458 = kfrho**2
785  t459 = t6*t458
786  t461 = t115*kfrho
787  t462 = t197*t461
788  t463 = t6*kfrhorho
789  t465 = t197*t463/0.2e1_dp
790  t466 = t99*t246
791  srhorho = t457*t459 + t462 - t465 + t466
792  s1rho = srho
793  t469 = mu**2
794  t470 = 0.1e1_dp/t202/t105*t469
795  t471 = t470*t101
796  t472 = srho*t103
797  t476 = s1rho*srho
798  fxrhorho = -0.8e1_dp*t471*t472*s1rho + 0.2e1_dp*t204* &
799  t476 + 0.2e1_dp*t204*s*srhorho
800  fx1rho = 0.2e1_dp*t204*s*s1rho
801  t487 = my_rho*ex_unifrhorho
802  t491 = my_rho*ex_unif1rho
803 
804  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
805  e_rho_rho(ii) = e_rho_rho(ii) + &
806  scale_ex*(ex_unif1rho*fx + ex_unif*fx1rho + &
807  ex_unifrho*fx + t487*fx + t208*fx1rho + ex_unif*fxrho + t491 &
808  *fxrho + t108*fxrhorho) + scale_ec*(e_c_u_01rho + t451*t191 &
809  + epsilon_cggarho + my_rho*epsilon_cggarhorho)
810  END IF
811 
812  t496 = t149*t6
813  t498 = t74*t115
814  tnorm_drhorho = -t496*k_srho/0.2e1_dp - t498/0.2e1_dp
815  t500 = t78*trho
816  t503 = t377*tnorm_drho
817  t506 = tnorm_drho*t186
818  t507 = t178*t506
819  t510 = t163*tnorm_drhorho
820  t513 = t91*trho
821  t517 = arho*tnorm_drho
822  t521 = a*tnorm_drhorho
823  t525 = t177*t186
824  t529 = t226*trho
825  t530 = t178*t529
826  t535 = t226*t186
827  t541 = a*trho
828  t548 = tnorm_drho*trho
829  t553 = 0.2e1_dp*t382*tnorm_drho + 0.2e1_dp*t541*tnorm_drho &
830  + 0.2e1_dp*t168*tnorm_drhorho + 0.8e1_dp*t417*t517 + &
831  0.12e2_dp*t426*t548 + 0.4e1_dp*t183*tnorm_drhorho
832  t556 = 0.2e1_dp*t500*t214 + 0.2e1_dp*t162*t503 - 0.2e1_dp &
833  *t162*t507 + 0.2e1_dp*t162*t510 + 0.6e1_dp*t175*t218* &
834  t513 + 0.2e1_dp*t217*t517*t91 + 0.2e1_dp*t217*t521*t91 - &
835  0.2e1_dp*t217*t218*t525 - 0.2e1_dp*t162*t530 - t175* &
836  t397*t226 + 0.2e1_dp*t175*t409*t535 - t175*t178*t553
837  t557 = gamma_var*t556
838  t559 = t439*t189
839  hnorm_drhorho = t557*t191 - t230*t559
840  t562 = t196*t6
841  snorm_drhorho = -t562*kfrho/0.2e1_dp - t98*t115/0.2e1_dp
842  t566 = snorm_drho*t103
843  fxnorm_drhorho = -0.8e1_dp*t471*t566*srho + 0.2e1_dp*t204 &
844  *srho*snorm_drho + 0.2e1_dp*t204*s*snorm_drhorho
845 
846  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
847  e_ndrho_rho(ii) = e_ndrho_rho(ii) + &
848  scale_ex*(ex_unif*fxnorm_drho + t208* &
849  fxnorm_drho + t108*fxnorm_drhorho) + scale_ec*(hnorm_drho + my_rho &
850  *hnorm_drhorho)
851  END IF
852 
853  t581 = tnorm_drho**2
854  t586 = a*t581
855  t590 = tnorm_drho*t226
856  t591 = t178*t590
857  t594 = t177*t226
858  t598 = t226**2
859  t605 = 0.2e1_dp*t586 + 0.12e2_dp*t426*t581
860  t609 = gamma_var*(0.2e1_dp*t78*t581*t85*t91 + 0.10e2_dp &
861  *t175*t586*t91 - 0.4e1_dp*t162*t591 - 0.4e1_dp*t217* &
862  t218*t594 + 0.2e1_dp*t175*t409*t598 - t175*t178*t605)
863  t611 = t229**2
864  t612 = gamma_var*t611
865  hnorm_drhonorm_drho = t609*t191 - t612*t439
866  t614 = snorm_drho**2
867  fxnorm_drhonorm_drho = -0.8e1_dp*t470*t101*t614*t103 + &
868  0.2e1_dp*t204*t614
869 
870  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
871  e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
872  scale_ex*t108*fxnorm_drhonorm_drho + &
873  scale_ec*my_rho*hnorm_drhonorm_drho
874  END IF
875 
876  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
877  rsrhorhorho = -0.5e1_dp/0.54e2_dp*t4/t111/t238/ &
878  t115*t627/t240/t114 + t4*t239/t240/my_rho/0.3e1_dp &
879  - t4*t113*t241/0.2e1_dp
880  rs2rho = rsrho
881  t645 = alpha_1_1*rsrhorho
882  t654 = t127*rs2rho/0.2e1_dp + beta_2_1*rs2rho + 0.3e1_dp/ &
883  0.2e1_dp*t131*rs2rho + t21*t19*rs2rho*t135
884  t656 = t124*t654*t139
885  t661 = t140*t654
886  t664 = rsrho*rs2rho
887  t669 = t21*t278
888  t670 = rs2rho*t281
889  t671 = t670*rsrho
890  t673 = t21*t19
891  t675 = -t266*t664/0.4e1_dp + t271 + t272 + 0.3e1_dp/0.4e1_dp &
892  *t273*t664 + t277 + t669*t671 + t286 - t673*t671
893  t685 = alpha_1_1*rs2rho
894  t693 = t675*t139
895  t701 = t297*t13
896  t702 = t701*t654
897  t714 = t267*rs2rho
898  t717 = rsrhorho*rsrho
899  t720 = rsrhorho*rs2rho
900  t740 = rs2rho/t280/rs*t267
901  t743 = rsrhorho*t281*rsrho
902  t748 = t670*rsrhorho
903  t758 = 0.3e1_dp/0.8e1_dp*beta_1_1/t14/t280*t714 - t266* &
904  t717/0.2e1_dp - t266*t720/0.4e1_dp + t127*rsrhorhorho/ &
905  0.2e1_dp + beta_2_1*rsrhorhorho - 0.3e1_dp/0.8e1_dp*beta_3_1* &
906  t265*t714 + 0.3e1_dp/0.2e1_dp*t273*t717 + 0.3e1_dp/ &
907  0.4e1_dp*t273*t720 + 0.3e1_dp/0.2e1_dp*t131*rsrhorhorho + &
908  t21*t278*t19*t740 + 0.2e1_dp*t669*t743 - 0.3e1_dp*t669* &
909  t740 + t669*t748 + t21*t19*rsrhorhorho*t135 - t673*t748 - &
910  0.2e1_dp*t673*t743 + 0.2e1_dp*t673*t740
911  t776 = a_1**2
912  e_c_u_0rhorhorho = -0.2e1_dp*t119*rsrhorhorho*t27 + t645* &
913  t656 + 0.2e1_dp*t645*t255 - 0.4e1_dp*t253*t259*t661 + &
914  0.2e1_dp*t253*t124*t675*t139 + t253*t294*t138*t297* &
915  t13*t654 - 0.2e1_dp*t685*t259*t261*t139 + 0.6e1_dp*t295 &
916  *t262*t654 - 0.4e1_dp*t260*t693*t138 - 0.3e1_dp*t11/ &
917  t293/t22*t261*t702 + t685*t124*t290*t139 - 0.2e1_dp* &
918  t260*t291*t654 + t125*t758*t139 + t295*t290*t702/ &
919  0.2e1_dp + t685*t294*t299/0.2e1_dp + t295*t675*t701*t138 &
920  + t11/t293/t123*t261/t296/t26/t776*t654/0.2e1_dp
921  e_c_u_0rho1rho = -t252 + t253*t656 + t685*t255 - 0.2e1_dp* &
922  t260*t661 + t125*t693 + t295*t138*t702/0.2e1_dp
923  e_c_u_01rhorho = e_c_u_0rho1rho
924  e_c_u_02rho = -0.2e1_dp*t119*rs2rho*t27 + t125*t654*t139
925  k_frhorhorho = 0.10e2_dp/0.27e2_dp*t1/t142/t114*t69
926  k_f2rho = kfrho
927  t801 = k_f**2
928  t809 = t309*k_frhorho
929  t812 = t238*k_f2rho
930  k_srho1rho = -t309*k_frho*t812/0.2e1_dp + t315
931  k_s1rhorho = k_srho1rho
932  k_s2rho = t146*k_f2rho*t5
933  t821 = t148**2
934  t825 = k_srho*k_s1rho
935  t831 = t6*k_srho1rho
936  trhorhorho = -0.3e1_dp*my_norm_drho/t821*t6*t825*k_s2rho - &
937  t318*t321*k_s1rho + t318*t831*k_s1rho + t318*t151* &
938  k_s1rhorho - t318*t321*k_s2rho - t150*t246*k_srho + t150* &
939  t115*k_srho1rho/0.2e1_dp + t318*t324*k_s2rho + t150*t115* &
940  k_srhorho/0.2e1_dp - t150*t6*(0.3e1_dp/0.4e1_dp/t73/ &
941  t801/t238*t310*t627*k_f2rho - t809*t238*k_frho - t809* &
942  t812/0.2e1_dp + t146*k_frhorhorho*t5)/0.2e1_dp - t318*t327 &
943  *k_s2rho - t150*t246*k_s1rho + t150*t115*k_s1rhorho/ &
944  0.2e1_dp - t150*t246*k_s2rho - 0.3e1_dp*t75*t241
945  t868 = t150*t115*k_s2rho/0.2e1_dp
946  trho1rho = t318*t151*k_s2rho + t323 - t150*t831/0.2e1_dp + &
947  t868 + t330
948  t1rhorho = t318*t331*k_s2rho + t329 - t150*t6*k_s1rhorho/ &
949  0.2e1_dp + t868 + t330
950  t2rho = -t150*t6*k_s2rho/0.2e1_dp - t153/0.2e1_dp
951  t877 = t155**2
952  t879 = beta/t877
953  t880 = t158**2
954  t885 = e_c_u_01rho*e_c_u_02rho
955  arhorhorho = 0.6e1_dp*t879/t880*e_c_u_0rho*t340*t80* &
956  t885 + 0.2e1_dp*t339*e_c_u_0rho1rho*t340*e_c_u_01rho - &
957  0.6e1_dp*t879*t338*t341*t885 + 0.2e1_dp*t339*t341* &
958  e_c_u_01rhorho + 0.2e1_dp*t339*e_c_u_0rhorho*t340* &
959  e_c_u_02rho + t157*t159*e_c_u_0rhorhorho*t80 - t348* &
960  e_c_u_0rhorho*e_c_u_02rho*t80 - t348*e_c_u_0rho1rho* &
961  e_c_u_01rho*t80 - t348*e_c_u_0rho*e_c_u_01rhorho*t80 + t879 &
962  *t159*t349*e_c_u_02rho*t80
963  arho1rho = 0.2e1_dp*t339*t341*e_c_u_02rho + t157*t159* &
964  e_c_u_0rho1rho*t80 - t348*e_c_u_0rho*e_c_u_02rho*t80
965  a1rhorho = 0.2e1_dp*t339*e_c_u_01rho*t340*e_c_u_02rho + &
966  t157*t159*e_c_u_01rhorho*t80 - t348*t885*t80
967  a2rho = t157*t159*e_c_u_02rho*t80
968  t940 = arho1rho*t83
969  t942 = 0.2e1_dp*t382*t2rho
970  t943 = a2rho*t
971  t945 = 0.2e1_dp*t943*trho
972  t946 = a*t2rho
973  t948 = 0.2e1_dp*t946*trho
974  t950 = 0.2e1_dp*t168*trho1rho
975  t951 = a2rho*t88
976  t954 = arho*t2rho
977  t967 = t940 + t942 + t945 + t948 + t950 + 0.2e1_dp*t951*arho + &
978  0.8e1_dp*t417*t954 + 0.2e1_dp*t179*arho1rho + 0.8e1_dp* &
979  t417*trho*a2rho + 0.12e2_dp*t426*trho*t2rho + 0.4e1_dp* &
980  t183*trho1rho
981  t976 = t78*t2rho
982  t980 = t78*t*t85
983  t982 = a2rho*t83
984  t984 = 0.2e1_dp*t168*t2rho
985  t989 = t982 + t984 + 0.2e1_dp*t179*a2rho + 0.4e1_dp*t183*t2rho
986  t990 = t369*t989
987  t994 = t982 + t984
988  t995 = t994*t177
989  t998 = arhorhorho*t83
990  t999 = arhorho*t
991  t1001 = 0.2e1_dp*t999*t2rho
992  t1004 = 0.2e1_dp*arho1rho*t*t1rho
993  t1005 = t954*t1rho
994  t1006 = 0.2e1_dp*t1005
995  t1008 = 0.2e1_dp*t382*t1rhorho
996  t1011 = 0.2e1_dp*a1rhorho*t*trho
997  t1012 = a1rho*t2rho
998  t1013 = t1012*trho
999  t1014 = 0.2e1_dp*t1013
1000  t1016 = 0.2e1_dp*t385*trho1rho
1001  t1017 = a2rho*t1rho
1002  t1018 = t1017*trho
1003  t1019 = 0.2e1_dp*t1018
1004  t1022 = 0.2e1_dp*a*t1rhorho*trho
1005  t1024 = 0.2e1_dp*t388*trho1rho
1006  t1026 = 0.2e1_dp*t943*trhorho
1007  t1028 = 0.2e1_dp*t946*trhorho
1008  t1030 = 0.2e1_dp*t168*trhorhorho
1009  t1031 = t998 + t1001 + t1004 + t1006 + t1008 + t1011 + t1014 + &
1010  t1016 + t1019 + t1022 + t1024 + t1026 + t1028 + t1030
1011  t1035 = t940 + t942 + t945 + t948 + t950
1012  t1042 = t369*t2rho
1013  t1046 = a1rhorho*t83
1014  t1048 = 0.2e1_dp*t385*t2rho
1015  t1050 = 0.2e1_dp*t943*t1rho
1016  t1052 = 0.2e1_dp*t946*t1rho
1017  t1054 = 0.2e1_dp*t168*t1rhorho
1018  t1055 = t1046 + t1048 + t1050 + t1052 + t1054
1019  t1060 = t78*t86
1020  t1061 = t176**2
1021  t1062 = 0.1e1_dp/t1061
1022  t1067 = -0.2e1_dp*t162*t178*t967*t1rho + 0.2e1_dp*t175* &
1023  t409*t967*t369 + 0.2e1_dp*t976*t374 + 0.4e1_dp*t980* &
1024  t408*trho*t990 - t175*t995*t432 + t78*t83*t1031*t91 - &
1025  t175*t1035*t177*t369 - 0.2e1_dp*t162*t995*t370 - &
1026  0.2e1_dp*t162*t397*t1042 + 0.2e1_dp*t162*t1055*t91* &
1027  trho - 0.6e1_dp*t1060*t1062*t186*t990
1028  t1093 = t1046 + t1048 + t1050 + t1052 + t1054 + 0.2e1_dp*t951* &
1029  a1rho + 0.8e1_dp*t417*t1012 + 0.2e1_dp*t179*a1rhorho + &
1030  0.8e1_dp*t417*t1017 + 0.12e2_dp*t426*t1rho*t2rho + &
1031  0.4e1_dp*t183*t1rhorho
1032  t1097 = t1rho*t989
1033  t1103 = trho*t989
1034  t1104 = t178*t1103
1035  t1106 = t994*t91
1036  t1109 = t171*t408
1037  t1115 = t976*t362 + t162*t361*trho1rho - t162*t178*t432 &
1038  *t2rho + t175*t994*t408*t410 - t162*t178*trho1rho*t369 &
1039  - t162*t178*trho*t1093 - t162*t397*t1097 + t175*t409* &
1040  t186*t1093 - t354*t1104 + t162*t1106*trhorho + t175*t1109 &
1041  *t990 + t175*t409*t432*t989
1042  t1118 = t1106*trho
1043  t1121 = t360*t408
1044  t1122 = t186*t989
1045  t1126 = t393*t177
1046  t1129 = t408*t186
1047  t1148 = t186*t2rho
1048  t1152 = 0.2e1_dp*t354*t1118 + 0.2e1_dp*t175*t1121*t1122 &
1049  - t175*t1126*t989 + 0.4e1_dp*t980*t1129*t1042 + 0.2e1_dp* &
1050  t976*t378 - t175*t404*t967 + 0.2e1_dp*t162*t377* &
1051  t1rhorho - 0.2e1_dp*t976*t401 + 0.2e1_dp*t78*t1rhorho*t164 &
1052  - 0.2e1_dp*t162*t404*t1103 - 0.2e1_dp*t162*t404*t1148
1053  t1157 = t393*t91
1054  t1167 = a2rho*t182
1055  t1187 = a1rho*t182
1056  t1196 = t87*t
1057  t1203 = t1008 + 0.8e1_dp*t417*trhorho*a2rho + 0.12e2_dp* &
1058  t426*trhorho*t2rho + 0.8e1_dp*t1167*t423 + 0.8e1_dp*t417* &
1059  trho*a1rhorho + 0.8e1_dp*t417*arho*t1rhorho + 0.8e1_dp* &
1060  t1167*t418 + 0.8e1_dp*t417*arho1rho*t1rho + 0.12e2_dp*t426 &
1061  *trho1rho*t1rho + 0.12e2_dp*t426*trho*t1rhorho + t998 + &
1062  0.8e1_dp*t1187*t954 + 0.8e1_dp*t417*trho1rho*a1rho + &
1063  0.8e1_dp*t417*arhorho*t2rho + 0.24e2_dp*t1196*t427*t2rho &
1064  + 0.2e1_dp*a1rhorho*t88*arho + t1014
1065  t1218 = t1016 + t1001 + 0.2e1_dp*t951*arhorho + t1026 + t1028 &
1066  + t1022 + t1011 + t1030 + 0.2e1_dp*t179*arhorhorho + 0.4e1_dp* &
1067  t183*trhorhorho + 0.2e1_dp*t414*arho1rho + t1004 + t1006 + &
1068  t1019 + t1024 + 0.24e2_dp*t84*t1018 + 0.24e2_dp*t84*t1005 + &
1069  0.24e2_dp*t84*t1013
1070  t1226 = t163*trho1rho
1071  t1249 = -0.2e1_dp*t162*t178*t186*t1rhorho + 0.2e1_dp* &
1072  t162*t1157*t2rho - t175*t178*(t1203 + t1218) + 0.2e1_dp* &
1073  t162*t1035*t91*t1rho + 0.2e1_dp*t354*t1226 - 0.2e1_dp* &
1074  t162*t995*t400 + 0.2e1_dp*t162*t163*trhorhorho - 0.2e1_dp &
1075  *t162*t178*trhorho*t989 - t175*t1055*t177*t186 - &
1076  0.2e1_dp*t976*t371 - t175*t397*t1093 + 0.4e1_dp*t980* &
1077  t1129*t1097
1078  t1262 = 0.2e1_dp*t162*t163*t2rho + t78*t83*t994*t91 - &
1079  t175*t178*t989
1080  t1263 = t439*t1262
1081  t1291 = 0.2e1_dp*t976*t164 + 0.2e1_dp*t162*t1118 - &
1082  0.2e1_dp*t162*t1104 + 0.2e1_dp*t162*t1226 + 0.2e1_dp*t162 &
1083  *t377*t2rho + t78*t83*t1035*t91 - t175*t397*t989 - &
1084  0.2e1_dp*t162*t178*t1148 - t175*t995*t186 + 0.2e1_dp* &
1085  t175*t409*t1122 - t175*t178*t967
1086  t1292 = gamma_var*t1291
1087  t1295 = 0.1e1_dp/t438/t94
1088  t1329 = 0.2e1_dp*t976*t440 + 0.2e1_dp*t162*t1106*t1rho - &
1089  0.2e1_dp*t162*t178*t1097 + 0.2e1_dp*t162*t163*t1rhorho &
1090  + 0.2e1_dp*t162*t361*t2rho + t78*t83*t1055*t91 - t175* &
1091  t404*t989 - 0.2e1_dp*t162*t178*t1042 - t175*t995*t369 + &
1092  0.2e1_dp*t175*t409*t990 - t175*t178*t1093
1093  kfrhorhorho = k_frhorhorho
1094  kf2rho = k_f2rho
1095  ex_unifrho1rho = ex_unifrhorho
1096  ex_unif1rhorho = ex_unifrho1rho
1097  ex_unif2rho = -0.3e1_dp/0.4e1_dp*t5*kf2rho
1098  t1342 = t195**2
1099  srho1rho = t457*t198*kf2rho + t462/0.2e1_dp - t465 + t197* &
1100  t115*kf2rho/0.2e1_dp + t466
1101  s1rhorho = srho1rho
1102  s2rho = -t197*t6*kf2rho/0.2e1_dp - t200/0.2e1_dp
1103  t1380 = t202**2
1104  t1385 = 0.1e1_dp/t1380*t469*mu*t101*s
1105  t1386 = kappa**2
1106  t1387 = 0.1e1_dp/t1386
1107  t1389 = s1rho*s2rho
1108  t1393 = t470*s
1109  fxrho1rho = -0.8e1_dp*t471*t472*s2rho + 0.2e1_dp*t204* &
1110  s2rho*srho + 0.2e1_dp*t204*s*srho1rho
1111  fx1rhorho = -0.8e1_dp*t471*s1rho*t103*s2rho + 0.2e1_dp* &
1112  t204*t1389 + 0.2e1_dp*t204*s*s1rhorho
1113  fx2rho = 0.2e1_dp*t204*s*s2rho
1114  ex_ldarhorhorho = ex_unif1rhorho*fx + ex_unif1rho*fx2rho + &
1115  ex_unif2rho*fx1rho + ex_unif*fx1rhorho + ex_unifrho1rho*fx + &
1116  ex_unifrho*fx2rho + ex_unifrhorho*fx - 0.3e1_dp/0.4e1_dp*my_rho &
1117  *t5*kfrhorhorho*fx + t487*fx2rho + ex_unifrho*fx1rho + my_rho &
1118  *ex_unifrho1rho*fx1rho + t208*fx1rhorho + ex_unif2rho*fxrho &
1119  + ex_unif*fxrho1rho + ex_unif1rho*fxrho + my_rho*ex_unif1rhorho* &
1120  fxrho + t491*fxrho1rho + ex_unif*fxrhorho + my_rho*ex_unif2rho* &
1121  fxrhorho + t108*(0.48e2_dp*t1385*srho*t1387*t1389 - &
1122  0.24e2_dp*t1393*t472*t1389 - 0.8e1_dp*t471*srho1rho*t103 &
1123  *s1rho - 0.8e1_dp*t471*t472*s1rhorho + 0.2e1_dp*t204* &
1124  s1rhorho*srho + 0.2e1_dp*t204*s1rho*srho1rho - 0.8e1_dp* &
1125  t471*srhorho*t103*s2rho + 0.2e1_dp*t204*s2rho*srhorho + &
1126  0.2e1_dp*t204*s*(-0.3e1_dp*my_norm_drho/t1342*t459*kf2rho &
1127  - t457*t115*t458 + 0.2e1_dp*t457*t463*kfrho - 0.2e1_dp* &
1128  t457*t461*kf2rho - 0.2e1_dp*t197*t246*kfrho + 0.3e1_dp/ &
1129  0.2e1_dp*t197*t115*kfrhorho + t457*t463*kf2rho - t197*t6 &
1130  *kfrhorhorho/0.2e1_dp - t197*t246*kf2rho - 0.3e1_dp*t99* &
1131  t241))
1132 
1133  e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + &
1134  scale_ex*ex_ldarhorhorho + scale_ec*( &
1135  e_c_u_01rhorho + gamma_var*t1329*t191 - t451*t1263 + &
1136  e_c_u_0rho1rho + t1292*t191 - t190*t1263 + epsilon_cggarhorho + &
1137  my_rho*(e_c_u_0rhorhorho + gamma_var*(t1067 + 0.2e1_dp*t1115 + &
1138  t1152 + t1249)*t191 - t436*t1263 - t1292*t449 + 0.2e1_dp* &
1139  t190*t1295*t448*t1262 - t190*t439*t1329))
1140 
1141  t1468 = t149*t115
1142  tnorm_drhorhorho = t317*t6*t825 + t1468*k_srho/0.2e1_dp - &
1143  t496*k_srhorho/0.2e1_dp + t1468*k_s1rho/0.2e1_dp + t74* &
1144  t246
1145  tnorm_drho1rho = -t496*k_s1rho/0.2e1_dp - t498/0.2e1_dp
1146  t1482 = a1rho*tnorm_drhorho
1147  t1492 = t78*t84
1148  t1493 = tnorm_drho*t177
1149  t1504 = t408*tnorm_drho
1150  t1505 = t1504*t410
1151  t1511 = a1rho*tnorm_drho
1152  t1515 = t226*t1rho
1153  t1521 = a*tnorm_drho1rho
1154  t1525 = 0.2e1_dp*t175*t409*t553*t369 + 0.2e1_dp*t217* &
1155  t1482*t91 - 0.2e1_dp*t162*t178*tnorm_drhorho*t369 + &
1156  0.2e1_dp*t354*t510 - 0.6e1_dp*t1492*t1493*t400 + 0.6e1_dp &
1157  *t175*t218*t91*trhorho + 0.2e1_dp*t162*t1157*tnorm_drho &
1158  + 0.4e1_dp*t980*t1505 - 0.6e1_dp*t1492*t1493*t370 + &
1159  0.6e1_dp*t175*t1511*t513 - 0.2e1_dp*t162*t397*t1515 - &
1160  0.2e1_dp*t354*t507 + 0.6e1_dp*t175*t1521*t513
1161  t1528 = arhorho*tnorm_drho
1162  t1532 = t177*t369
1163  t1545 = t78*t417
1164  t1565 = 0.2e1_dp*t385*tnorm_drho + 0.2e1_dp*t388* &
1165  tnorm_drho + 0.2e1_dp*t168*tnorm_drho1rho + 0.8e1_dp*t417* &
1166  t1511 + 0.12e2_dp*t426*tnorm_drho*t1rho + 0.4e1_dp*t183* &
1167  tnorm_drho1rho
1168  t1573 = t91*t1rho
1169  t1584 = -0.2e1_dp*t354*t530 + 0.2e1_dp*t217*t1528*t91 - &
1170  0.2e1_dp*t217*t517*t1532 - 0.2e1_dp*t217*t1511*t525 + &
1171  0.2e1_dp*t162*t377*tnorm_drho1rho + 0.2e1_dp*t162*t361* &
1172  tnorm_drhorho + 0.4e1_dp*t1545*t1505 + 0.2e1_dp*t217*a* &
1173  tnorm_drhorhorho*t91 + 0.2e1_dp*t175*t409*t1565*t186 - &
1174  0.2e1_dp*t217*t521*t1532 + 0.6e1_dp*t175*t521*t1573 - &
1175  0.6e1_dp*t1060*t1062*t226*t410 + 0.6e1_dp*t175*t517* &
1176  t1573
1177  t1608 = t408*t226
1178  t1612 = t163*tnorm_drho1rho
1179  t1628 = -0.2e1_dp*t162*t404*t506 + 0.2e1_dp*t354*t503 - &
1180  0.2e1_dp*t162*t178*tnorm_drho*t432 - 0.2e1_dp*t162*t178 &
1181  *t553*t1rho - 0.2e1_dp*t162*t404*t529 - 0.2e1_dp*t162* &
1182  t178*t1565*trho - t175*t1126*t226 + 0.4e1_dp*t980*t1608 &
1183  *t370 + 0.2e1_dp*t500*t1612 + 0.12e2_dp*t78*t168* &
1184  tnorm_drho*t91*t427 + 0.2e1_dp*t162*t163*tnorm_drhorhorho &
1185  - t175*t404*t553 + 0.2e1_dp*t175*t1121*t535
1186  t1629 = arho*tnorm_drho1rho
1187  t1633 = tnorm_drho*t369
1188  t1646 = t361*tnorm_drho
1189  t1652 = t226*t369
1190  t1660 = t178*t1633
1191  t1672 = t418*tnorm_drho
1192  t1676 = t423*tnorm_drho
1193  t1715 = 0.2e1_dp*t999*tnorm_drho + 0.2e1_dp*t1672 + 0.2e1_dp &
1194  *t382*tnorm_drho1rho + 0.2e1_dp*t1676 + 0.2e1_dp*a*trhorho &
1195  *tnorm_drho + 0.2e1_dp*t541*tnorm_drho1rho + 0.2e1_dp*t385* &
1196  tnorm_drhorho + 0.2e1_dp*t388*tnorm_drhorho + 0.2e1_dp*t168* &
1197  tnorm_drhorhorho + 0.8e1_dp*t1187*t517 + 0.24e2_dp*t84* &
1198  t1672 + 0.8e1_dp*t417*t1629 + 0.8e1_dp*t417*t1528 + &
1199  0.24e2_dp*t84*t1676 + 0.24e2_dp*t1196*t548*t1rho + &
1200  0.12e2_dp*t426*tnorm_drho1rho*trho + 0.12e2_dp*t426* &
1201  tnorm_drho*trhorho + 0.8e1_dp*t417*t1482 + 0.12e2_dp*t426* &
1202  tnorm_drhorho*t1rho + 0.4e1_dp*t183*tnorm_drhorhorho
1203  t1722 = 0.2e1_dp*t217*t1629*t91 - 0.2e1_dp*t162*t397* &
1204  t1633 - t175*t397*t1565 - 0.2e1_dp*t162*t178*t226* &
1205  trhorho + 0.2e1_dp*t78*trhorho*t214 + 0.2e1_dp*t500*t1646 &
1206  - 0.2e1_dp*t217*t1521*t525 + 0.2e1_dp*t175*t1109*t1652 - &
1207  0.2e1_dp*t162*t178*tnorm_drho1rho*t186 - 0.2e1_dp*t500* &
1208  t1660 + 0.4e1_dp*t980*t1608*t400 + 0.2e1_dp*t175*t409* &
1209  t226*t432 - t175*t178*t1715 - 0.2e1_dp*t217*t218*t177* &
1210  t432
1211  t1758 = 0.2e1_dp*t354*t214 + 0.2e1_dp*t162*t1646 - &
1212  0.2e1_dp*t162*t1660 + 0.2e1_dp*t162*t1612 + 0.6e1_dp*t175 &
1213  *t218*t1573 + 0.2e1_dp*t217*t1511*t91 + 0.2e1_dp*t217* &
1214  t1521*t91 - 0.2e1_dp*t217*t218*t1532 - 0.2e1_dp*t162* &
1215  t178*t1515 - t175*t404*t226 + 0.2e1_dp*t175*t409*t1652 - &
1216  t175*t178*t1565
1217  t1759 = gamma_var*t1758
1218  t1761 = t1295*t189
1219  snorm_drho1rho = snorm_drhorho
1220  t1797 = snorm_drhorho*t103
1221  fxnorm_drho1rho = -0.8e1_dp*t471*t566*s1rho + 0.2e1_dp* &
1222  t204*s1rho*snorm_drho + 0.2e1_dp*t204*s*snorm_drho1rho
1223 
1224  e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) + &
1225  scale_ex*(ex_unif1rho*fxnorm_drho + &
1226  ex_unif*fxnorm_drho1rho + ex_unifrho*fxnorm_drho + t487* &
1227  fxnorm_drho + t208*fxnorm_drho1rho + ex_unif*fxnorm_drhorho + &
1228  t491*fxnorm_drhorho + t108*(0.48e2_dp*t1385*snorm_drho* &
1229  t1387*t476 - 0.24e2_dp*t1393*t566*t476 - 0.8e1_dp*t471* &
1230  snorm_drho1rho*t103*srho - 0.8e1_dp*t471*t566*srhorho + &
1231  0.2e1_dp*t204*srhorho*snorm_drho + 0.2e1_dp*t204*srho* &
1232  snorm_drho1rho - 0.8e1_dp*t471*t1797*s1rho + 0.2e1_dp*t204* &
1233  s1rho*snorm_drhorho + 0.2e1_dp*t204*s*(t456*t6*t458 + &
1234  t196*t115*kfrho - t562*kfrhorho/0.2e1_dp + t98*t246))) + &
1235  scale_ec*(t1759*t191 - t230*t449 + hnorm_drhorho + my_rho*( &
1236  gamma_var*(t1525 + t1584 + t1628 + t1722)*t191 - t557*t449 - &
1237  t1759*t559 + 0.2e1_dp*t230*t1761*t448 - t230*t439*t435))
1238 
1239  t1838 = t1504*t535
1240  t1851 = arho*t581
1241  t1878 = 0.4e1_dp*t175*t409*t226*t553 - 0.2e1_dp*t162* &
1242  t178*t605*trho - 0.4e1_dp*t217*t218*t177*t553 + 0.8e1_dp &
1243  *t1545*t1838 + 0.2e1_dp*t78*t581*t171*t91 - 0.4e1_dp* &
1244  t162*t397*t590 - 0.10e2_dp*t175*t586*t525 - t175*t178* &
1245  (0.2e1_dp*t1851 + 0.4e1_dp*t521*tnorm_drho + 0.24e2_dp*t84* &
1246  t1851 + 0.24e2_dp*t1196*t581*trho + 0.24e2_dp*t426* &
1247  tnorm_drhorho*tnorm_drho) - 0.12e2_dp*t1492*t1493*t529 + &
1248  0.8e1_dp*t980*t1838 - 0.4e1_dp*t162*t178*tnorm_drhorho* &
1249  t226 + 0.20e2_dp*t162*t586*t513
1250  t1889 = t78*t581
1251  t1907 = t85*t1062
1252  t1922 = -0.4e1_dp*t500*t591 + 0.2e1_dp*t175*t409*t605* &
1253  t186 + 0.4e1_dp*t162*t409*t598*trho - 0.2e1_dp*t1889* &
1254  t187 + 0.2e1_dp*t175*t1109*t598 + 0.20e2_dp*t175*t218* &
1255  t91*tnorm_drhorho + 0.10e2_dp*t175*t1851*t91 - 0.4e1_dp* &
1256  t217*t517*t594 - t175*t397*t605 - 0.6e1_dp*t175*t1907* &
1257  t598*t186 - 0.4e1_dp*t162*t178*tnorm_drho*t553 + 0.4e1_dp &
1258  *t78*tnorm_drhorho*t214 - 0.4e1_dp*t217*t521*t594
1259  t1927 = t439*t229
1260  t1933 = t614*t1387
1261  t1937 = t614*t103
1262 
1263  e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) + &
1264  scale_ex*(ex_unif* &
1265  fxnorm_drhonorm_drho + t208*fxnorm_drhonorm_drho + t108*( &
1266  0.48e2_dp*t1385*t1933*srho - 0.24e2_dp*t1393*t1937*srho &
1267  - 0.16e2_dp*t471*t1797*snorm_drho + 0.4e1_dp*t204* &
1268  snorm_drhorho*snorm_drho)) + scale_ec*(hnorm_drhonorm_drho + my_rho &
1269  *(gamma_var*(t1878 + t1922)*t191 - t609*t559 - 0.2e1_dp* &
1270  t557*t1927 + 0.2e1_dp*t612*t1761))
1271 
1272  t2norm_drho = tnorm_drho
1273  t1952 = t226*t2norm_drho
1274  t1964 = 0.2e1_dp*t168*t2norm_drho + 0.4e1_dp*t183*t2norm_drho
1275  t1965 = t178*t1964
1276  t1968 = t226*t1964
1277  t1969 = t1504*t1968
1278  t1972 = a*t2norm_drho
1279  t1990 = 0.2e1_dp*t1972*tnorm_drho + 0.12e2_dp*t426* &
1280  tnorm_drho*t2norm_drho
1281  t2020 = t177*t1964
1282  t2024 = t91*t2norm_drho
1283  t2028 = t78*t2norm_drho
1284  t2031 = -0.20e2_dp*t1492*t1493*t1952 + 0.4e1_dp*t162* &
1285  t409*t598*t2norm_drho - 0.2e1_dp*t1889*t1965 + 0.8e1_dp* &
1286  t1545*t1969 + 0.4e1_dp*t217*t1972*t408*t598 - 0.6e1_dp* &
1287  t175*t1907*t598*t1964 - 0.2e1_dp*t217*t1972*t177*t605 &
1288  - 0.4e1_dp*t217*t218*t177*t1990 + 0.4e1_dp*t175*t409* &
1289  t1990*t226 - 0.24e2_dp*t78*t182*t85*t177*t87*t581* &
1290  t2norm_drho + 0.2e1_dp*t175*t409*t605*t1964 + 0.8e1_dp* &
1291  t980*t1969 - 0.2e1_dp*t162*t178*t605*t2norm_drho - &
1292  0.4e1_dp*t162*t178*t1990*tnorm_drho - 0.10e2_dp*t175* &
1293  t586*t2020 + 0.24e2_dp*t162*t586*t2024 - 0.4e1_dp*t2028* &
1294  t591
1295  t2041 = 0.2e1_dp*t162*t163*t2norm_drho + 0.2e1_dp*t217* &
1296  t1972*t91 - t175*t1965
1297  s2norm_drho = snorm_drho
1298 
1299  e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) + &
1300  scale_ex*t108*(0.48e2_dp* &
1301  t1385*t1933*s2norm_drho - 0.24e2_dp*t1393*t1937* &
1302  s2norm_drho) + scale_ec*my_rho*(gamma_var*t2031*t191 - t609* &
1303  t439*t2041 - 0.2e1_dp*gamma_var*(0.2e1_dp*t2028*t214 + &
1304  0.10e2_dp*t175*t218*t2024 - 0.2e1_dp*t162*t178* &
1305  tnorm_drho*t1964 - 0.2e1_dp*t217*t218*t2020 - 0.2e1_dp* &
1306  t162*t178*t1952 - 0.2e1_dp*t217*t1972*t594 + 0.2e1_dp* &
1307  t175*t409*t1968 - t175*t178*t1990)*t1927 + 0.2e1_dp*t612 &
1308  *t1295*t2041)
1309  END IF
1310  END IF
1311  END DO
1312 !$OMP END DO
1313  END SELECT
1314 
1315  END SUBROUTINE pbe_lda_calc
1316 
1317 ! **************************************************************************************************
1318 !> \brief evaluates the becke 88 exchange functional for lsd
1319 !> \param rho_set ...
1320 !> \param deriv_set ...
1321 !> \param grad_deriv ...
1322 !> \param pbe_params ...
1323 !> \author fawzi
1324 ! **************************************************************************************************
1325  SUBROUTINE pbe_lsd_eval(rho_set, deriv_set, grad_deriv, pbe_params)
1326  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1327  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
1328  INTEGER, INTENT(in) :: grad_deriv
1329  TYPE(section_vals_type), POINTER :: pbe_params
1330 
1331  CHARACTER(len=*), PARAMETER :: routinen = 'pbe_lsd_eval'
1332 
1333  INTEGER :: handle, npoints, param
1334  INTEGER, DIMENSION(2, 3) :: bo
1335  REAL(kind=dp) :: epsilon_rho, scale_ec, scale_ex
1336  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, &
1337  e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndrb, e_ndrb_ndrb, e_ndrb_rb, e_ra, &
1338  e_ra_ra, e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, norm_drhob, rhoa, rhob
1339  TYPE(xc_derivative_type), POINTER :: deriv
1340 
1341  CALL timeset(routinen, handle)
1342  NULLIFY (deriv)
1343 
1344  CALL xc_rho_set_get(rho_set, &
1345  rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
1346  norm_drhob=norm_drhob, norm_drho=norm_drho, &
1347  rho_cutoff=epsilon_rho, &
1348  local_bounds=bo)
1349  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
1350 
1351  dummy => rhoa
1352  e_0 => dummy
1353  e_ra => dummy
1354  e_rb => dummy
1355  e_ndra_ra => dummy
1356  e_ndrb_rb => dummy
1357  e_ndr_ndr => dummy
1358  e_ndra_ndra => dummy
1359  e_ndrb_ndrb => dummy
1360  e_ndr => dummy
1361  e_ndra => dummy
1362  e_ndrb => dummy
1363  e_ra_ra => dummy
1364  e_ra_rb => dummy
1365  e_rb_rb => dummy
1366  e_ndr_ra => dummy
1367  e_ndr_rb => dummy
1368 
1369  IF (grad_deriv >= 0) THEN
1370  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
1371  allocate_deriv=.true.)
1372  CALL xc_derivative_get(deriv, deriv_data=e_0)
1373  END IF
1374  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1375  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
1376  allocate_deriv=.true.)
1377  CALL xc_derivative_get(deriv, deriv_data=e_ra)
1378  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
1379  allocate_deriv=.true.)
1380  CALL xc_derivative_get(deriv, deriv_data=e_rb)
1381  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1382  allocate_deriv=.true.)
1383  CALL xc_derivative_get(deriv, deriv_data=e_ndr)
1384  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
1385  allocate_deriv=.true.)
1386  CALL xc_derivative_get(deriv, deriv_data=e_ndra)
1387  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
1388  allocate_deriv=.true.)
1389  CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
1390  END IF
1391  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1392  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1393  allocate_deriv=.true.)
1394  CALL xc_derivative_get(deriv, deriv_data=e_ra_ra)
1395  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1396  allocate_deriv=.true.)
1397  CALL xc_derivative_get(deriv, deriv_data=e_ra_rb)
1398  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1399  allocate_deriv=.true.)
1400  CALL xc_derivative_get(deriv, deriv_data=e_rb_rb)
1401  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
1402  allocate_deriv=.true.)
1403  CALL xc_derivative_get(deriv, deriv_data=e_ndr_ra)
1404  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
1405  allocate_deriv=.true.)
1406  CALL xc_derivative_get(deriv, deriv_data=e_ndr_rb)
1407  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
1408  allocate_deriv=.true.)
1409  CALL xc_derivative_get(deriv, deriv_data=e_ndra_ra)
1410  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
1411  allocate_deriv=.true.)
1412  CALL xc_derivative_get(deriv, deriv_data=e_ndrb_rb)
1413  deriv => xc_dset_get_derivative(deriv_set, &
1414  [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
1415  CALL xc_derivative_get(deriv, deriv_data=e_ndr_ndr)
1416  deriv => xc_dset_get_derivative(deriv_set, &
1417  [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.true.)
1418  CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndra)
1419  deriv => xc_dset_get_derivative(deriv_set, &
1420  [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.true.)
1421  CALL xc_derivative_get(deriv, deriv_data=e_ndrb_ndrb)
1422  END IF
1423  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
1424  ! to do
1425  END IF
1426 
1427  CALL section_vals_val_get(pbe_params, "scale_c", r_val=scale_ec)
1428  CALL section_vals_val_get(pbe_params, "scale_x", r_val=scale_ex)
1429  CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
1430 
1431 !$OMP PARALLEL DEFAULT(NONE),&
1432 !$OMP SHARED(rhoa,rhob,norm_drho,norm_drhoa,norm_drhob,e_0,e_ra,e_rb,e_ndra_ra),&
1433 !$OMP SHARED(e_ndrb_rb,e_ndr_ndr,e_ndra_ndra,e_ndrb_ndrb,e_ndr,e_ndra,e_ndrb),&
1434 !$OMP SHARED(e_ra_ra,e_ra_rb,e_rb_rb,e_ndr_ra,e_ndr_rb,grad_deriv,npoints),&
1435 !$OMP SHARED(epsilon_rho,param,scale_ec,scale_ex)
1436  CALL pbe_lsd_calc( &
1437  rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa, &
1438  norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb, &
1439  e_ra_ndra=e_ndra_ra, &
1440  e_rb_ndrb=e_ndrb_rb, e_ndr_ndr=e_ndr_ndr, &
1441  e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, e_ndr=e_ndr, &
1442  e_ndra=e_ndra, e_ndrb=e_ndrb, e_ra_ra=e_ra_ra, &
1443  e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, e_ra_ndr=e_ndr_ra, &
1444  e_rb_ndr=e_ndr_rb, &
1445  grad_deriv=grad_deriv, npoints=npoints, &
1446  epsilon_rho=epsilon_rho, &
1447  param=param, scale_ec=scale_ec, scale_ex=scale_ex)
1448 !$OMP END PARALLEL
1449 
1450  CALL timestop(handle)
1451  END SUBROUTINE pbe_lsd_eval
1452 
1453 ! **************************************************************************************************
1454 !> \brief low level calculation of the pbe exchange-correlation functional for lsd
1455 !> \param rhoa ,rhob: alpha or beta spin density
1456 !> \param rhob ...
1457 !> \param norm_drho ...
1458 !> \param norm_drhoa ,norm_drhob,norm_drho: || grad rhoa |||,| grad rhoa ||,
1459 !> || grad (rhoa+rhob) ||
1460 !> \param norm_drhob ...
1461 !> \param e_0 adds to it the local value of the functional
1462 !> \param e_ra derivative of the functional wrt. to the variables
1463 !> named where the * is
1464 !> \param e_rb ...
1465 !> \param e_ra_ndra ...
1466 !> \param e_rb_ndrb ...
1467 !> \param e_ndr_ndr ...
1468 !> \param e_ndra_ndra ...
1469 !> \param e_ndrb_ndrb ...
1470 !> \param e_ndr ...
1471 !> \param e_ndra ...
1472 !> \param e_ndrb ...
1473 !> \param e_ra_ra ...
1474 !> \param e_ra_rb ...
1475 !> \param e_rb_rb ...
1476 !> \param e_ra_ndr ...
1477 !> \param e_rb_ndr ...
1478 !> \param grad_deriv ...
1479 !> \param npoints ...
1480 !> \param epsilon_rho ...
1481 !> \param param ...
1482 !> \param scale_ec ...
1483 !> \param scale_ex ...
1484 !> \author fawzi
1485 ! **************************************************************************************************
1486  SUBROUTINE pbe_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, &
1487  e_0, e_ra, e_rb, e_ra_ndra, e_rb_ndrb, e_ndr_ndr, &
1488  e_ndra_ndra, e_ndrb_ndrb, e_ndr, &
1489  e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ra_ndr, e_rb_ndr, &
1490  grad_deriv, npoints, epsilon_rho, param, scale_ec, scale_ex)
1491  REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, norm_drhoa, &
1492  norm_drhob
1493  REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ra_ndra, e_rb_ndrb, &
1494  e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb, e_ndr, e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, &
1495  e_ra_ndr, e_rb_ndr
1496  INTEGER, INTENT(in) :: grad_deriv, npoints
1497  REAL(kind=dp), INTENT(in) :: epsilon_rho
1498  INTEGER, INTENT(in) :: param
1499  REAL(kind=dp), INTENT(in) :: scale_ec, scale_ex
1500 
1501  INTEGER :: ii
1502  REAL(kind=dp) :: a, a1rhoa, a1rhob, a_1, a_2, a_3, alpha_1_1, alpha_1_2, alpha_1_3, alpha_c, &
1503  alpha_c1rhoa, alpha_c1rhob, alpha_crhoa, alpha_crhob, arhoa, arhoarhoa, arhoarhob, arhob, &
1504  arhobrhob, beta, beta_1_1, beta_1_2, beta_1_3, beta_2_1, beta_2_2, beta_2_3, beta_3_1, &
1505  beta_3_2, beta_3_3, beta_4_1, beta_4_2, beta_4_3, chi, chirhoa, chirhoarhoa, chirhoarhob, &
1506  chirhob, chirhobrhob, e_c_u_0, e_c_u_01rhoa, e_c_u_01rhob, e_c_u_0rhoa, e_c_u_0rhoarhoa, &
1507  e_c_u_0rhoarhob, e_c_u_0rhob, e_c_u_0rhobrhob, e_c_u_1rhoa, e_c_u_1rhob, epsilon_c_unif, &
1508  epsilon_c_unif1rhoa, epsilon_c_unif1rhob
1509  REAL(kind=dp) :: epsilon_c_unifrhoa, epsilon_c_unifrhoarhoa, epsilon_c_unifrhoarhob, &
1510  epsilon_c_unifrhob, epsilon_c_unifrhobrhob, epsilon_cgga, epsilon_cggarhoa, &
1511  epsilon_cggarhob, ex_unif_a, ex_unif_a1rhoa, ex_unif_arhoa, ex_unif_b, ex_unif_b1rhob, &
1512  ex_unif_brhob, f, f1rhoa, f1rhob, frhoa, frhoarhoa, frhoarhob, frhob, frhobrhob, fx_a, &
1513  fx_a1rhoa, fx_anorm_drhoa, fx_arhoa, fx_b, fx_b1rhob, fx_bnorm_drhob, fx_brhob, &
1514  gamma_var, hnorm_drho, k_frhoa, k_frhoarhoa, k_frhoarhob, k_frhob, k_s, k_s1rhoa, &
1515  k_s1rhob, k_srhoa, k_srhob, kappa, kf_a, kf_arhoa, kf_arhoarhoa, kf_b, kf_brhob, &
1516  kf_brhobrhob, mu
1517  REAL(kind=dp) :: my_norm_drho, my_norm_drhoa, my_norm_drhob, my_rho, my_rhoa, my_rhob, p_1, &
1518  p_2, p_3, phi, phi1rhoa, phi1rhob, phirhoa, phirhoarhoa, phirhoarhob, phirhob, &
1519  phirhobrhob, rs, rsrhoa, rsrhoarhoa, rsrhoarhob, rsrhob, rsrhobrhob, s_a, s_a1rhoa, &
1520  s_anorm_drhoa, s_arhoa, s_b, s_b1rhob, s_bnorm_drhob, s_brhob, t, t1, t100, t1000, t1001, &
1521  t101, t102, t103, t1031, t1033, t1038, t104, t105, t1050, t1069, t107, t1071, t108, t110, &
1522  t1104, t1106, t111, t112, t113, t115, t116, t1164, t118, t119, t1193, t12, t122, t1228, &
1523  t123, t1231, t124, t125, t126, t1269, t1283, t1285, t1286, t1288, t129
1524  REAL(kind=dp) :: t1291, t1293, t130, t1304, t131, t1327, t133, t1330, t1342, t1346, t135, &
1525  t137, t1377, t14, t140, t1411, t142, t143, t1440, t146, t1462, t1465, t147, t148, t1482, &
1526  t1489, t1492, t15, t150, t1501, t1514, t1520, t1529, t153, t1550, t1552, t1555, t156, &
1527  t1588, t1590, t1591, t1599, t1602, t1603, t162, t163, t1630, t1632, t1635, t1638, t164, &
1528  t1640, t165, t167, t1674, t1677, t1680, t1698, t171, t1711, t1712, t1713, t1739, t1741, &
1529  t1743, t1748, t175, t176, t1765, t1766, t1767, t177, t178, t179, t18, t1801, t1829, t183, &
1530  t1830, t1831, t1865, t187, t1871, t1876, t1888, t190, t1901, t191
1531  REAL(kind=dp) :: t192, t1922, t194, t1949, t198, t199, t1rhoa, t1rhob, t2, t20, t200, t201, &
1532  t205, t21, t211, t212, t213, t215, t219, t22, t220, t221, t222, t226, t23, t232, t233, &
1533  t234, t240, t242, t244, t245, t246, t248, t249, t250, t252, t254, t256, t257, t259, t262, &
1534  t266, t268, t269, t27, t272, t273, t274, t277, t278, t28, t280, t281, t282, t284, t285, &
1535  t286, t289, t293, t294, t297, t298, t299, t3, t302, t303, t305, t306, t310, t311, t312, &
1536  t313, t314, t317, t318, t32, t321, t324, t325, t326, t329, t335, t336, t337, t34, t340, &
1537  t341, t344, t346, t350, t368, t38, t382, t39, t396, t4, t40
1538  REAL(kind=dp) :: t403, t405, t407, t409, t41, t410, t411, t413, t415, t417, t427, t429, &
1539  t432, t436, t439, t442, t444, t445, t45, t453, t456, t457, t46, t460, t466, t467, t468, &
1540  t471, t472, t475, t477, t481, t488, t493, t494, t5, t50, t502, t505, t506, t518, t519, &
1541  t52, t523, t525, t536, t538, t543, t544, t548, t549, t550, t556, t56, t561, t563, t564, &
1542  t57, t576, t578, t579, t58, t580, t588, t59, t590, t595, t596, t6, t600, t606, t611, &
1543  t624, t626, t627, t628, t63, t636, t638, t64, t643, t644, t648, t654, t659, t66, t672, &
1544  t674, t675, t676, t681, t682, t687, t69, t7, t70, t705, t708, t71, t711
1545  REAL(kind=dp) :: t72, t726, t73, t733, t736, t74, t745, t75, t750, t755, t763, t767, t768, &
1546  t77, t775, t776, t779, t78, t785, t789, t79, t795, t798, t8, t80, t801, t812, t82, t820, &
1547  t821, t822, t823, t825, t828, t839, t84, t840, t85, t851, t858, t865, t867, t868, t87, &
1548  t876, t879, t88, t880, t9, t90, t904, t908, t909, t91, t911, t914, t917, t919, t92, t924, &
1549  t93, t936, t94, t944, t95, t953, t959, t96, t962, t965, t966, t967, t97, t98, t985, t998, &
1550  t999, tnorm_drho, trhoa, trhoanorm_drho, trhoarhoa, trhoarhob, trhob, trhobnorm_drho, &
1551  trhobrhob
1552 
1553 ! Parametrization of PBE
1554 
1555  SELECT CASE (param)
1556  ! Original PBE
1557  CASE (xc_pbe_orig)
1558  kappa = 0.804e0_dp
1559  beta = 0.66725e-1_dp
1560  mu = beta*pi**2/0.3e1_dp
1561  ! Revised PBE (revPBE)
1562  CASE (xc_pbe_rev)
1563  kappa = 0.1245e1_dp
1564  beta = 0.66725e-1_dp
1565  mu = beta*pi**2/0.3e1_dp
1566  ! PBE for solids (PBEsol)
1567  CASE (xc_pbe_sol)
1568  kappa = 0.804e0_dp
1569  beta = 0.46e-1_dp
1570  mu = 0.1e2_dp/0.81e2_dp
1571  CASE default
1572  cpabort("Unsupported parametrization")
1573  END SELECT
1574 
1575  gamma_var = (0.1e1_dp - log(0.2e1_dp))/pi**2
1576  p_1 = 0.10e1_dp
1577  a_1 = 0.31091e-1_dp
1578  alpha_1_1 = 0.21370e0_dp
1579  beta_1_1 = 0.75957e1_dp
1580  beta_2_1 = 0.35876e1_dp
1581  beta_3_1 = 0.16382e1_dp
1582  beta_4_1 = 0.49294e0_dp
1583  p_2 = 0.10e1_dp
1584  a_2 = 0.15545e-1_dp
1585  alpha_1_2 = 0.20548e0_dp
1586  beta_1_2 = 0.141189e2_dp
1587  beta_2_2 = 0.61977e1_dp
1588  beta_3_2 = 0.33662e1_dp
1589  beta_4_2 = 0.62517e0_dp
1590  p_3 = 0.10e1_dp
1591  a_3 = 0.16887e-1_dp
1592  alpha_1_3 = 0.11125e0_dp
1593  beta_1_3 = 0.10357e2_dp
1594  beta_2_3 = 0.36231e1_dp
1595  beta_3_3 = 0.88026e0_dp
1596  beta_4_3 = 0.49671e0_dp
1597  t3 = 3**(0.1e1_dp/0.3e1_dp)
1598  t4 = 4**(0.1e1_dp/0.3e1_dp)
1599  t5 = t4**2
1600  t6 = t3*t5
1601  t7 = 0.1e1_dp/pi
1602  t90 = pi**2
1603 
1604  SELECT CASE (grad_deriv)
1605  CASE default
1606 !$OMP DO
1607  DO ii = 1, npoints
1608  my_rhoa = max(rhoa(ii), 0.0_dp)
1609  my_rhob = max(rhob(ii), 0.0_dp)
1610  my_rho = my_rhoa + my_rhob
1611  IF (my_rho > epsilon_rho) THEN
1612  my_rhoa = max(epsilon(0.0_dp)*1.e4_dp, my_rhoa)
1613  my_rhob = max(epsilon(0.0_dp)*1.e4_dp, my_rhob)
1614  my_rho = my_rhoa + my_rhob
1615  my_norm_drho = norm_drho(ii)
1616  my_norm_drhoa = norm_drhoa(ii)
1617  my_norm_drhob = norm_drhob(ii)
1618 
1619  t1 = my_rhoa - my_rhob
1620  t2 = 0.1e1_dp/my_rho
1621  chi = t1*t2
1622  t8 = t7*t2
1623  t9 = t8**(0.1e1_dp/0.3e1_dp)
1624  rs = t6*t9/0.4e1_dp
1625  t12 = 0.1e1_dp + alpha_1_1*rs
1626  t14 = 0.1e1_dp/a_1
1627  t15 = sqrt(rs)
1628  t18 = t15*rs
1629  t20 = p_1 + 0.1e1_dp
1630  t21 = rs**t20
1631  t22 = beta_4_1*t21
1632  t23 = beta_1_1*t15 + beta_2_1*rs + beta_3_1*t18 + t22
1633  t27 = 0.1e1_dp + t14/t23/0.2e1_dp
1634  t28 = log(t27)
1635  e_c_u_0 = -0.2e1_dp*a_1*t12*t28
1636  t32 = 0.1e1_dp + alpha_1_2*rs
1637  t34 = 0.1e1_dp/a_2
1638  t38 = p_2 + 0.1e1_dp
1639  t39 = rs**t38
1640  t40 = beta_4_2*t39
1641  t41 = beta_1_2*t15 + beta_2_2*rs + beta_3_2*t18 + t40
1642  t45 = 0.1e1_dp + t34/t41/0.2e1_dp
1643  t46 = log(t45)
1644  t50 = 0.1e1_dp + alpha_1_3*rs
1645  t52 = 0.1e1_dp/a_3
1646  t56 = p_3 + 0.1e1_dp
1647  t57 = rs**t56
1648  t58 = beta_4_3*t57
1649  t59 = beta_1_3*t15 + beta_2_3*rs + beta_3_3*t18 + t58
1650  t63 = 0.1e1_dp + t52/t59/0.2e1_dp
1651  t64 = log(t63)
1652  alpha_c = 0.2e1_dp*a_3*t50*t64
1653  t66 = 2**(0.1e1_dp/0.3e1_dp)
1654  t69 = 1/(2*t66 - 2)
1655  t70 = 0.1e1_dp + chi
1656  t71 = t70**(0.1e1_dp/0.3e1_dp)
1657  t72 = t71*t70
1658  t73 = 0.1e1_dp - chi
1659  t74 = t73**(0.1e1_dp/0.3e1_dp)
1660  t75 = t74*t73
1661  f = (t72 + t75 - 0.2e1_dp)*t69
1662  t77 = alpha_c*f
1663  t78 = 0.9e1_dp/0.8e1_dp/t69
1664  t79 = chi**2
1665  t80 = t79**2
1666  t82 = t78*(0.1e1_dp - t80)
1667  t84 = -0.2e1_dp*a_2*t32*t46 - e_c_u_0
1668  t85 = t84*f
1669  epsilon_c_unif = e_c_u_0 + t77*t82 + t85*t80
1670  t87 = t71**2
1671  t88 = t74**2
1672  phi = t87/0.2e1_dp + t88/0.2e1_dp
1673  t91 = t90*my_rho
1674  t92 = t91**(0.1e1_dp/0.3e1_dp)
1675  t93 = t3*t92*t7
1676  t94 = sqrt(t93)
1677  k_s = 0.2e1_dp*t94
1678  t95 = 0.1e1_dp/phi
1679  t96 = my_norm_drho*t95
1680  t97 = 0.1e1_dp/k_s
1681  t98 = t97*t2
1682  t = t96*t98/0.2e1_dp
1683  t100 = 0.1e1_dp/gamma_var
1684  t101 = beta*t100
1685  t102 = epsilon_c_unif*t100
1686  t103 = phi**2
1687  t104 = t103*phi
1688  t105 = 0.1e1_dp/t104
1689  t107 = exp(-t102*t105)
1690  t108 = t107 - 0.1e1_dp
1691  a = t101/t108
1692  t110 = gamma_var*t104
1693  t111 = t**2
1694  t112 = a*t111
1695  t113 = 0.1e1_dp + t112
1696  t115 = a**2
1697  t116 = t111**2
1698  t118 = 0.1e1_dp + t112 + t115*t116
1699  t119 = 0.1e1_dp/t118
1700  t122 = 0.1e1_dp + t101*t111*t113*t119
1701  t123 = log(t122)
1702  epsilon_cgga = epsilon_c_unif + t110*t123
1703  t124 = t3*t66
1704  t125 = t90*my_rhoa
1705  t126 = t125**(0.1e1_dp/0.3e1_dp)
1706  kf_a = t124*t126
1707  ex_unif_a = -0.3e1_dp/0.4e1_dp*t7*kf_a
1708  t129 = 0.1e1_dp/kf_a
1709  t130 = my_norm_drhoa*t129
1710  t131 = 0.1e1_dp/my_rhoa
1711  s_a = t130*t131/0.2e1_dp
1712  t133 = s_a**2
1713  t135 = 0.1e1_dp/kappa
1714  t137 = 0.1e1_dp + mu*t133*t135
1715  fx_a = 0.1e1_dp + kappa - kappa/t137
1716  t140 = my_rhoa*ex_unif_a
1717  t142 = t90*my_rhob
1718  t143 = t142**(0.1e1_dp/0.3e1_dp)
1719  kf_b = t124*t143
1720  ex_unif_b = -0.3e1_dp/0.4e1_dp*t7*kf_b
1721  t146 = 0.1e1_dp/kf_b
1722  t147 = my_norm_drhob*t146
1723  t148 = 0.1e1_dp/my_rhob
1724  s_b = t147*t148/0.2e1_dp
1725  t150 = s_b**2
1726  t153 = 0.1e1_dp + mu*t150*t135
1727  fx_b = 0.1e1_dp + kappa - kappa/t153
1728  t156 = my_rhob*ex_unif_b
1729 
1730  IF (grad_deriv >= 0) THEN
1731  e_0(ii) = e_0(ii) + &
1732  scale_ex*(0.2e1_dp*t140*fx_a + 0.2e1_dp*t156*fx_b) &
1733  /0.2e1_dp + scale_ec*my_rho*epsilon_cgga
1734  END IF
1735 
1736  t162 = my_rho**2
1737  t163 = 0.1e1_dp/t162
1738  t164 = t1*t163
1739  chirhoa = t2 - t164
1740  t165 = t9**2
1741  t167 = 0.1e1_dp/t165*t7
1742  rsrhoa = -t6*t167*t163/0.12e2_dp
1743  t171 = a_1*alpha_1_1
1744  t175 = t23**2
1745  t176 = 0.1e1_dp/t175
1746  t177 = t12*t176
1747  t178 = 0.1e1_dp/t15
1748  t179 = beta_1_1*t178
1749  t183 = beta_3_1*t15
1750  t187 = 0.1e1_dp/rs
1751  t190 = t179*rsrhoa/0.2e1_dp + beta_2_1*rsrhoa + 0.3e1_dp/ &
1752  0.2e1_dp*t183*rsrhoa + t22*t20*rsrhoa*t187
1753  t191 = 0.1e1_dp/t27
1754  t192 = t190*t191
1755  e_c_u_0rhoa = -0.2e1_dp*t171*rsrhoa*t28 + t177*t192
1756  t194 = a_2*alpha_1_2
1757  t198 = t41**2
1758  t199 = 0.1e1_dp/t198
1759  t200 = t32*t199
1760  t201 = beta_1_2*t178
1761  t205 = beta_3_2*t15
1762  t211 = t201*rsrhoa/0.2e1_dp + beta_2_2*rsrhoa + 0.3e1_dp/ &
1763  0.2e1_dp*t205*rsrhoa + t40*t38*rsrhoa*t187
1764  t212 = 0.1e1_dp/t45
1765  t213 = t211*t212
1766  e_c_u_1rhoa = -0.2e1_dp*t194*rsrhoa*t46 + t200*t213
1767  t215 = a_3*alpha_1_3
1768  t219 = t59**2
1769  t220 = 0.1e1_dp/t219
1770  t221 = t50*t220
1771  t222 = beta_1_3*t178
1772  t226 = beta_3_3*t15
1773  t232 = t222*rsrhoa/0.2e1_dp + beta_2_3*rsrhoa + 0.3e1_dp/ &
1774  0.2e1_dp*t226*rsrhoa + t58*t56*rsrhoa*t187
1775  t233 = 0.1e1_dp/t63
1776  t234 = t232*t233
1777  alpha_crhoa = 0.2e1_dp*t215*rsrhoa*t64 - t221*t234
1778  frhoa = (0.4e1_dp/0.3e1_dp*t71*chirhoa - 0.4e1_dp/0.3e1_dp &
1779  *t74*chirhoa)*t69
1780  t240 = alpha_crhoa*f
1781  t242 = alpha_c*frhoa
1782  t244 = t79*chi
1783  t245 = t78*t244
1784  t246 = t245*chirhoa
1785  t248 = 0.4e1_dp*t77*t246
1786  t249 = e_c_u_1rhoa - e_c_u_0rhoa
1787  t250 = t249*f
1788  t252 = t84*frhoa
1789  t254 = t244*chirhoa
1790  t256 = 0.4e1_dp*t85*t254
1791  epsilon_c_unifrhoa = e_c_u_0rhoa + t240*t82 + t242*t82 - t248 &
1792  + t250*t80 + t252*t80 + t256
1793  t257 = 0.1e1_dp/t71
1794  t259 = 0.1e1_dp/t74
1795  phirhoa = t257*chirhoa/0.3e1_dp - t259*chirhoa/0.3e1_dp
1796  t262 = t92**2
1797  k_frhoa = t3/t262*t90/0.3e1_dp
1798  t266 = 0.1e1_dp/t94
1799  k_srhoa = t266*k_frhoa*t7
1800  t268 = 0.1e1_dp/t103
1801  t269 = my_norm_drho*t268
1802  t272 = k_s**2
1803  t273 = 0.1e1_dp/t272
1804  t274 = t273*t2
1805  t277 = t97*t163
1806  t278 = t96*t277
1807  trhoa = -t269*t98*phirhoa/0.2e1_dp - t96*t274*k_srhoa/ &
1808  0.2e1_dp - t278/0.2e1_dp
1809  t280 = t108**2
1810  t281 = 0.1e1_dp/t280
1811  t282 = epsilon_c_unifrhoa*t100
1812  t284 = t103**2
1813  t285 = 0.1e1_dp/t284
1814  t286 = t285*phirhoa
1815  t289 = -t282*t105 + 0.3e1_dp*t102*t286
1816  arhoa = -t101*t281*t289*t107
1817  t293 = gamma_var*t103
1818  t294 = t123*phirhoa
1819  t297 = t101*t
1820  t298 = t113*t119
1821  t299 = t298*trhoa
1822  t302 = arhoa*t111
1823  t303 = a*t
1824  t305 = 0.2e1_dp*t303*trhoa
1825  t306 = t302 + t305
1826  t310 = t101*t111
1827  t311 = t118**2
1828  t312 = 0.1e1_dp/t311
1829  t313 = t113*t312
1830  t314 = a*t116
1831  t317 = t111*t
1832  t318 = t115*t317
1833  t321 = t302 + t305 + 0.2e1_dp*t314*arhoa + 0.4e1_dp*t318*trhoa
1834  t324 = 0.2e1_dp*t297*t299 + t101*t111*t306*t119 - t310* &
1835  t313*t321
1836  t325 = 0.1e1_dp/t122
1837  t326 = t324*t325
1838  epsilon_cggarhoa = epsilon_c_unifrhoa + 0.3e1_dp*t293*t294 + &
1839  t110*t326
1840  t329 = t126**2
1841  kf_arhoa = t124/t329*t90/0.3e1_dp
1842  ex_unif_arhoa = -0.3e1_dp/0.4e1_dp*t7*kf_arhoa
1843  t335 = kf_a**2
1844  t336 = 0.1e1_dp/t335
1845  t337 = my_norm_drhoa*t336
1846  t340 = my_rhoa**2
1847  t341 = 0.1e1_dp/t340
1848  s_arhoa = -t337*t131*kf_arhoa/0.2e1_dp - t130*t341/0.2e1_dp
1849  t344 = t137**2
1850  t346 = 0.1e1_dp/t344*mu
1851  fx_arhoa = 0.2e1_dp*t346*s_a*s_arhoa
1852  t350 = my_rhoa*ex_unif_arhoa
1853 
1854  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1855  e_ra(ii) = e_ra(ii) + &
1856  scale_ex*(0.2e1_dp*ex_unif_a*fx_a + 0.2e1_dp* &
1857  t350*fx_a + 0.2e1_dp*t140*fx_arhoa)/0.2e1_dp + scale_ec*( &
1858  epsilon_cgga + my_rho*epsilon_cggarhoa)
1859  END IF
1860 
1861  chirhob = -t2 - t164
1862  rsrhob = rsrhoa
1863  t368 = t179*rsrhob/0.2e1_dp + beta_2_1*rsrhob + 0.3e1_dp/ &
1864  0.2e1_dp*t183*rsrhob + t22*t20*rsrhob*t187
1865  e_c_u_0rhob = -0.2e1_dp*t171*rsrhob*t28 + t177*t368*t191
1866  t382 = t201*rsrhob/0.2e1_dp + beta_2_2*rsrhob + 0.3e1_dp/ &
1867  0.2e1_dp*t205*rsrhob + t40*t38*rsrhob*t187
1868  e_c_u_1rhob = -0.2e1_dp*t194*rsrhob*t46 + t200*t382*t212
1869  t396 = t222*rsrhob/0.2e1_dp + beta_2_3*rsrhob + 0.3e1_dp/ &
1870  0.2e1_dp*t226*rsrhob + t58*t56*rsrhob*t187
1871  alpha_crhob = 0.2e1_dp*t215*rsrhob*t64 - t221*t396*t233
1872  frhob = (0.4e1_dp/0.3e1_dp*t71*chirhob - 0.4e1_dp/0.3e1_dp &
1873  *t74*chirhob)*t69
1874  t403 = alpha_crhob*f
1875  t405 = alpha_c*frhob
1876  t407 = t245*chirhob
1877  t409 = 0.4e1_dp*t77*t407
1878  t410 = e_c_u_1rhob - e_c_u_0rhob
1879  t411 = t410*f
1880  t413 = t84*frhob
1881  t415 = t244*chirhob
1882  t417 = 0.4e1_dp*t85*t415
1883  epsilon_c_unifrhob = e_c_u_0rhob + t403*t82 + t405*t82 - t409 &
1884  + t411*t80 + t413*t80 + t417
1885  phirhob = t257*chirhob/0.3e1_dp - t259*chirhob/0.3e1_dp
1886  k_frhob = k_frhoa
1887  k_srhob = t266*k_frhob*t7
1888  trhob = -t269*t98*phirhob/0.2e1_dp - t96*t274*k_srhob/ &
1889  0.2e1_dp - t278/0.2e1_dp
1890  t427 = epsilon_c_unifrhob*t100
1891  t429 = t285*phirhob
1892  t432 = -t427*t105 + 0.3e1_dp*t102*t429
1893  arhob = -t101*t281*t432*t107
1894  t436 = t123*phirhob
1895  t439 = t298*trhob
1896  t442 = arhob*t111
1897  t444 = 0.2e1_dp*t303*trhob
1898  t445 = t442 + t444
1899  t453 = t442 + t444 + 0.2e1_dp*t314*arhob + 0.4e1_dp*t318*trhob
1900  t456 = 0.2e1_dp*t297*t439 + t101*t111*t445*t119 - t310* &
1901  t313*t453
1902  t457 = t456*t325
1903  epsilon_cggarhob = epsilon_c_unifrhob + 0.3e1_dp*t293*t436 + &
1904  t110*t457
1905  t460 = t143**2
1906  kf_brhob = t124/t460*t90/0.3e1_dp
1907  ex_unif_brhob = -0.3e1_dp/0.4e1_dp*t7*kf_brhob
1908  t466 = kf_b**2
1909  t467 = 0.1e1_dp/t466
1910  t468 = my_norm_drhob*t467
1911  t471 = my_rhob**2
1912  t472 = 0.1e1_dp/t471
1913  s_brhob = -t468*t148*kf_brhob/0.2e1_dp - t147*t472/0.2e1_dp
1914  t475 = t153**2
1915  t477 = 0.1e1_dp/t475*mu
1916  fx_brhob = 0.2e1_dp*t477*s_b*s_brhob
1917  t481 = my_rhob*ex_unif_brhob
1918 
1919  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1920  e_rb(ii) = e_rb(ii) + &
1921  scale_ex*(0.2e1_dp*ex_unif_b*fx_b + 0.2e1_dp* &
1922  t481*fx_b + 0.2e1_dp*t156*fx_brhob)/0.2e1_dp + scale_ec*( &
1923  epsilon_cgga + my_rho*epsilon_cggarhob)
1924  END IF
1925 
1926  t488 = t95*t97
1927  tnorm_drho = t488*t2/0.2e1_dp
1928  t493 = t101*t317
1929  t494 = a*tnorm_drho
1930  t502 = 0.2e1_dp*t303*tnorm_drho + 0.4e1_dp*t318*tnorm_drho
1931  t505 = 0.2e1_dp*t297*t298*tnorm_drho + 0.2e1_dp*t493* &
1932  t494*t119 - t310*t313*t502
1933  t506 = t505*t325
1934  hnorm_drho = t110*t506
1935 
1936  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1937  e_ndr(ii) = e_ndr(ii) + &
1938  scale_ec*my_rho*hnorm_drho
1939  END IF
1940 
1941  s_anorm_drhoa = t129*t131/0.2e1_dp
1942  fx_anorm_drhoa = 0.2e1_dp*t346*s_a*s_anorm_drhoa
1943 
1944  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1945  e_ndra(ii) = e_ndra(ii) + &
1946  scale_ex*t140*fx_anorm_drhoa
1947  END IF
1948 
1949  s_bnorm_drhob = t146*t148/0.2e1_dp
1950  fx_bnorm_drhob = 0.2e1_dp*t477*s_b*s_bnorm_drhob
1951 
1952  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1953  e_ndrb(ii) = e_ndrb(ii) + &
1954  scale_ex*t156*fx_bnorm_drhob
1955  END IF
1956 
1957  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1958  t518 = 0.1e1_dp/t162/my_rho
1959  t519 = t1*t518
1960  chirhoarhoa = -0.2e1_dp*t163 + 0.2e1_dp*t519
1961  t523 = 0.1e1_dp/t90
1962  t525 = t162**2
1963  rsrhoarhoa = -t6/t165/t8*t523/t525/0.18e2_dp + &
1964  t6*t167*t518/0.6e1_dp
1965  t536 = alpha_1_1*rsrhoa
1966  t538 = t176*t190*t191
1967  t543 = t12/t175/t23
1968  t544 = t190**2
1969  t548 = 0.1e1_dp/t18
1970  t549 = beta_1_1*t548
1971  t550 = rsrhoa**2
1972  t556 = beta_3_1*t178
1973  t561 = t20**2
1974  t563 = rs**2
1975  t564 = 0.1e1_dp/t563
1976  t576 = t175**2
1977  t578 = t12/t576
1978  t579 = t27**2
1979  t580 = 0.1e1_dp/t579
1980  e_c_u_0rhoarhoa = -0.2e1_dp*t171*rsrhoarhoa*t28 + 0.2e1_dp* &
1981  t536*t538 - 0.2e1_dp*t543*t544*t191 + t177*(-t549*t550 &
1982  /0.4e1_dp + t179*rsrhoarhoa/0.2e1_dp + beta_2_1*rsrhoarhoa + &
1983  0.3e1_dp/0.4e1_dp*t556*t550 + 0.3e1_dp/0.2e1_dp*t183* &
1984  rsrhoarhoa + t22*t561*t550*t564 + t22*t20*rsrhoarhoa* &
1985  t187 - t22*t20*t550*t564)*t191 + t578*t544*t580*t14/ &
1986  0.2e1_dp
1987  e_c_u_01rhoa = e_c_u_0rhoa
1988  t588 = alpha_1_2*rsrhoa
1989  t590 = t199*t211*t212
1990  t595 = t32/t198/t41
1991  t596 = t211**2
1992  t600 = beta_1_2*t548
1993  t606 = beta_3_2*t178
1994  t611 = t38**2
1995  t624 = t198**2
1996  t626 = t32/t624
1997  t627 = t45**2
1998  t628 = 0.1e1_dp/t627
1999  t636 = alpha_1_3*rsrhoa
2000  t638 = t220*t232*t233
2001  t643 = t50/t219/t59
2002  t644 = t232**2
2003  t648 = beta_1_3*t548
2004  t654 = beta_3_3*t178
2005  t659 = t56**2
2006  t672 = t219**2
2007  t674 = t50/t672
2008  t675 = t63**2
2009  t676 = 0.1e1_dp/t675
2010  alpha_c1rhoa = alpha_crhoa
2011  t681 = 0.1e1_dp/t87
2012  t682 = chirhoa**2
2013  t687 = 0.1e1_dp/t88
2014  frhoarhoa = (0.4e1_dp/0.9e1_dp*t681*t682 + 0.4e1_dp/ &
2015  0.3e1_dp*t71*chirhoarhoa + 0.4e1_dp/0.9e1_dp*t687*t682 - &
2016  0.4e1_dp/0.3e1_dp*t74*chirhoarhoa)*t69
2017  f1rhoa = frhoa
2018  t705 = alpha_c1rhoa*f
2019  t708 = alpha_c*f1rhoa
2020  t711 = t78*t79
2021  t726 = e_c_u_1rhoa - e_c_u_01rhoa
2022  t733 = t726*f
2023  t736 = t84*f1rhoa
2024  t745 = -0.4e1_dp*t77*t245*chirhoarhoa + (-0.2e1_dp*t194* &
2025  rsrhoarhoa*t46 + 0.2e1_dp*t588*t590 - 0.2e1_dp*t595*t596* &
2026  t212 + t200*(-t600*t550/0.4e1_dp + t201*rsrhoarhoa/ &
2027  0.2e1_dp + beta_2_2*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t606*t550 &
2028  + 0.3e1_dp/0.2e1_dp*t205*rsrhoarhoa + t40*t611*t550* &
2029  t564 + t40*t38*rsrhoarhoa*t187 - t40*t38*t550*t564)* &
2030  t212 + t626*t596*t628*t34/0.2e1_dp - e_c_u_0rhoarhoa)*f* &
2031  t80 + t249*f1rhoa*t80 + 0.4e1_dp*t250*t254 + t726*frhoa* &
2032  t80 + t84*frhoarhoa*t80 + 0.4e1_dp*t252*t254 + 0.4e1_dp* &
2033  t733*t254 + 0.4e1_dp*t736*t254 + 0.12e2_dp*t85*t79*t682 &
2034  + 0.4e1_dp*t85*t244*chirhoarhoa
2035  epsilon_c_unifrhoarhoa = e_c_u_0rhoarhoa + (0.2e1_dp*t215* &
2036  rsrhoarhoa*t64 - 0.2e1_dp*t636*t638 + 0.2e1_dp*t643*t644* &
2037  t233 - t221*(-t648*t550/0.4e1_dp + t222*rsrhoarhoa/ &
2038  0.2e1_dp + beta_2_3*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t654*t550 &
2039  + 0.3e1_dp/0.2e1_dp*t226*rsrhoarhoa + t58*t659*t550* &
2040  t564 + t58*t56*rsrhoarhoa*t187 - t58*t56*t550*t564)* &
2041  t233 - t674*t644*t676*t52/0.2e1_dp)*f*t82 + alpha_crhoa &
2042  *f1rhoa*t82 - 0.4e1_dp*t240*t246 + alpha_c1rhoa*frhoa*t82 &
2043  + alpha_c*frhoarhoa*t82 - 0.4e1_dp*t242*t246 - 0.4e1_dp* &
2044  t705*t246 - 0.4e1_dp*t708*t246 - 0.12e2_dp*t77*t711*t682 &
2045  + t745
2046  epsilon_c_unif1rhoa = e_c_u_01rhoa + t705*t82 + t708*t82 - &
2047  t248 + t733*t80 + t736*t80 + t256
2048  t750 = 0.1e1_dp/t72
2049  t755 = 0.1e1_dp/t75
2050  phirhoarhoa = -t750*t682/0.9e1_dp + t257*chirhoarhoa/ &
2051  0.3e1_dp - t755*t682/0.9e1_dp - t259*chirhoarhoa/0.3e1_dp
2052  phi1rhoa = phirhoa
2053  t763 = t90**2
2054  k_frhoarhoa = -0.2e1_dp/0.9e1_dp*t3/t262/t91*t763
2055  t767 = 0.1e1_dp/t94/t93
2056  t768 = k_frhoa**2
2057  k_s1rhoa = k_srhoa
2058  t775 = my_norm_drho*t105*t97
2059  t776 = t2*phirhoa
2060  t779 = t269*t273
2061  t785 = t269*t277*phirhoa/0.2e1_dp
2062  t789 = t2*k_srhoa
2063  t795 = t96/t272/k_s
2064  t798 = t273*t163
2065  t801 = t96*t798*k_srhoa/0.2e1_dp
2066  t812 = t96*t97*t518
2067  trhoarhoa = t775*t776*phi1rhoa + t779*t776*k_s1rhoa/ &
2068  0.2e1_dp + t785 - t269*t98*phirhoarhoa/0.2e1_dp + t779*t789 &
2069  *phi1rhoa/0.2e1_dp + t795*t789*k_s1rhoa + t801 - t96*t274* &
2070  (-t767*t768*t523/0.2e1_dp + t266*k_frhoarhoa*t7)/ &
2071  0.2e1_dp + t269*t277*phi1rhoa/0.2e1_dp + t96*t798*k_s1rhoa &
2072  /0.2e1_dp + t812
2073  t1rhoa = -t269*t98*phi1rhoa/0.2e1_dp - t96*t274*k_s1rhoa &
2074  /0.2e1_dp - t278/0.2e1_dp
2075  t820 = t101/t280/t108
2076  t821 = t107**2
2077  t822 = t289*t821
2078  t823 = epsilon_c_unif1rhoa*t100
2079  t825 = t285*phi1rhoa
2080  t828 = -t823*t105 + 0.3e1_dp*t102*t825
2081  t839 = 0.1e1_dp/t284/phi
2082  t840 = t839*phirhoa
2083  t851 = t101*t281
2084  arhoarhoa = 0.2e1_dp*t820*t822*t828 - t101*t281*( &
2085  -epsilon_c_unifrhoarhoa*t100*t105 + 0.3e1_dp*t282*t825 + &
2086  0.3e1_dp*t823*t286 - 0.12e2_dp*t102*t840*phi1rhoa + &
2087  0.3e1_dp*t102*t285*phirhoarhoa)*t107 - t851*t289*t828* &
2088  t107
2089  a1rhoa = -t101*t281*t828*t107
2090  t858 = gamma_var*phi
2091  t865 = a1rhoa*t111
2092  t867 = 0.2e1_dp*t303*t1rhoa
2093  t868 = t865 + t867
2094  t876 = t865 + t867 + 0.2e1_dp*t314*a1rhoa + 0.4e1_dp*t318*t1rhoa
2095  t879 = 0.2e1_dp*t297*t298*t1rhoa + t101*t111*t868*t119 &
2096  - t310*t313*t876
2097  t880 = t879*t325
2098  t904 = t306*t119
2099  t908 = arhoarhoa*t111
2100  t909 = arhoa*t
2101  t911 = 0.2e1_dp*t909*t1rhoa
2102  t914 = 0.2e1_dp*a1rhoa*t*trhoa
2103  t917 = 0.2e1_dp*a*t1rhoa*trhoa
2104  t919 = 0.2e1_dp*t303*trhoarhoa
2105  t924 = t306*t312
2106  t936 = t113/t311/t118
2107  t944 = a*t317
2108  t953 = t115*t111
2109  t959 = t908 + t911 + t914 + t917 + t919 + 0.2e1_dp*a1rhoa*t116 &
2110  *arhoa + 0.8e1_dp*t944*arhoa*t1rhoa + 0.2e1_dp*t314* &
2111  arhoarhoa + 0.8e1_dp*t944*trhoa*a1rhoa + 0.12e2_dp*t953* &
2112  trhoa*t1rhoa + 0.4e1_dp*t318*trhoarhoa
2113  t962 = 0.2e1_dp*t101*t1rhoa*t299 + 0.2e1_dp*t297*t868* &
2114  t119*trhoa - 0.2e1_dp*t297*t313*trhoa*t876 + 0.2e1_dp* &
2115  t297*t298*trhoarhoa + 0.2e1_dp*t297*t904*t1rhoa + t101* &
2116  t111*(t908 + t911 + t914 + t917 + t919)*t119 - t310*t924* &
2117  t876 - 0.2e1_dp*t297*t313*t321*t1rhoa - t310*t868*t312* &
2118  t321 + 0.2e1_dp*t310*t936*t321*t876 - t310*t313*t959
2119  t965 = t122**2
2120  t966 = 0.1e1_dp/t965
2121  t967 = t324*t966
2122  kf_arhoarhoa = -0.2e1_dp/0.9e1_dp*t124/t329/t125*t763
2123  ex_unif_a1rhoa = ex_unif_arhoa
2124  t985 = kf_arhoa**2
2125  s_a1rhoa = s_arhoa
2126  t998 = mu**2
2127  t999 = 0.1e1_dp/t344/t137*t998
2128  t1000 = t999*t133
2129  t1001 = s_arhoa*t135
2130  fx_a1rhoa = 0.2e1_dp*t346*s_a*s_a1rhoa
2131 
2132  e_ra_ra(ii) = e_ra_ra(ii) + &
2133  scale_ex*(0.2e1_dp*ex_unif_a1rhoa*fx_a + &
2134  0.2e1_dp*ex_unif_a*fx_a1rhoa + 0.2e1_dp*ex_unif_arhoa*fx_a - &
2135  0.3e1_dp/0.2e1_dp*my_rhoa*t7*kf_arhoarhoa*fx_a + 0.2e1_dp* &
2136  t350*fx_a1rhoa + 0.2e1_dp*ex_unif_a*fx_arhoa + 0.2e1_dp*my_rhoa &
2137  *ex_unif_a1rhoa*fx_arhoa + 0.2e1_dp*t140*(-0.8e1_dp*t1000 &
2138  *t1001*s_a1rhoa + 0.2e1_dp*t346*s_a1rhoa*s_arhoa + 0.2e1_dp &
2139  *t346*s_a*(my_norm_drhoa/t335/kf_a*t131*t985 + t337* &
2140  t341*kf_arhoa - t337*t131*kf_arhoarhoa/0.2e1_dp + t130/ &
2141  t340/my_rhoa)))/0.2e1_dp + scale_ec*(epsilon_c_unif1rhoa + &
2142  0.3e1_dp*t293*t123*phi1rhoa + t110*t880 + epsilon_cggarhoa + &
2143  my_rho*(epsilon_c_unifrhoarhoa + 0.6e1_dp*t858*t294*phi1rhoa + &
2144  0.3e1_dp*t293*t880*phirhoa + 0.3e1_dp*t293*t123* &
2145  phirhoarhoa + 0.3e1_dp*t293*t326*phi1rhoa + t110*t962*t325 &
2146  - t110*t967*t879))
2147 
2148  chirhoarhob = 0.2e1_dp*t519
2149  rsrhoarhob = rsrhoarhoa
2150  t1031 = t176*t368*t191
2151  t1033 = alpha_1_1*rsrhob
2152  t1038 = rsrhoa*rsrhob
2153  t1050 = rsrhob*t564*rsrhoa
2154  e_c_u_0rhoarhob = -0.2e1_dp*t171*rsrhoarhob*t28 + t536* &
2155  t1031 + t1033*t538 - 0.2e1_dp*t543*t192*t368 + t177*(-t549 &
2156  *t1038/0.4e1_dp + t179*rsrhoarhob/0.2e1_dp + beta_2_1* &
2157  rsrhoarhob + 0.3e1_dp/0.4e1_dp*t556*t1038 + 0.3e1_dp/ &
2158  0.2e1_dp*t183*rsrhoarhob + t22*t561*t1050 + t22*t20* &
2159  rsrhoarhob*t187 - t22*t20*t1050)*t191 + t578*t190*t580* &
2160  t14*t368/0.2e1_dp
2161  t1069 = t199*t382*t212
2162  t1071 = alpha_1_2*rsrhob
2163  t1104 = t220*t396*t233
2164  t1106 = alpha_1_3*rsrhob
2165  frhoarhob = (0.4e1_dp/0.9e1_dp*t681*chirhoa*chirhob + &
2166  0.4e1_dp/0.3e1_dp*t71*chirhoarhob + 0.4e1_dp/0.9e1_dp*t687 &
2167  *chirhoa*chirhob - 0.4e1_dp/0.3e1_dp*t74*chirhoarhob)* &
2168  t69
2169  t1164 = t79*chirhoa*chirhob
2170  t1193 = -0.4e1_dp*t77*t245*chirhoarhob + (-0.2e1_dp*t194* &
2171  rsrhoarhob*t46 + t588*t1069 + t1071*t590 - 0.2e1_dp*t595* &
2172  t213*t382 + t200*(-t600*t1038/0.4e1_dp + t201*rsrhoarhob/ &
2173  0.2e1_dp + beta_2_2*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t606* &
2174  t1038 + 0.3e1_dp/0.2e1_dp*t205*rsrhoarhob + t40*t611*t1050 &
2175  + t40*t38*rsrhoarhob*t187 - t40*t38*t1050)*t212 + t626 &
2176  *t211*t628*t34*t382/0.2e1_dp - e_c_u_0rhoarhob)*f*t80 + &
2177  t249*frhob*t80 + 0.4e1_dp*t250*t415 + t410*frhoa*t80 + &
2178  t84*frhoarhob*t80 + 0.4e1_dp*t252*t415 + 0.4e1_dp*t411* &
2179  t254 + 0.4e1_dp*t413*t254 + 0.12e2_dp*t85*t1164 + 0.4e1_dp* &
2180  t85*t244*chirhoarhob
2181  epsilon_c_unifrhoarhob = e_c_u_0rhoarhob + (0.2e1_dp*t215* &
2182  rsrhoarhob*t64 - t636*t1104 - t1106*t638 + 0.2e1_dp*t643* &
2183  t234*t396 - t221*(-t648*t1038/0.4e1_dp + t222*rsrhoarhob/ &
2184  0.2e1_dp + beta_2_3*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t654* &
2185  t1038 + 0.3e1_dp/0.2e1_dp*t226*rsrhoarhob + t58*t659*t1050 &
2186  + t58*t56*rsrhoarhob*t187 - t58*t56*t1050)*t233 - t674 &
2187  *t232*t676*t52*t396/0.2e1_dp)*f*t82 + alpha_crhoa* &
2188  frhob*t82 - 0.4e1_dp*t240*t407 + alpha_crhob*frhoa*t82 + &
2189  alpha_c*frhoarhob*t82 - 0.4e1_dp*t242*t407 - 0.4e1_dp*t403 &
2190  *t246 - 0.4e1_dp*t405*t246 - 0.12e2_dp*t77*t78*t1164 + &
2191  t1193
2192  phirhoarhob = -t750*chirhoa*chirhob/0.9e1_dp + t257* &
2193  chirhoarhob/0.3e1_dp - t755*chirhoa*chirhob/0.9e1_dp - t259 &
2194  *chirhoarhob/0.3e1_dp
2195  k_frhoarhob = k_frhoarhoa
2196  t1228 = t269*t277*phirhob/0.2e1_dp
2197  t1231 = t96*t798*k_srhob/0.2e1_dp
2198  trhoarhob = t775*t776*phirhob + t779*t776*k_srhob/ &
2199  0.2e1_dp + t785 - t269*t98*phirhoarhob/0.2e1_dp + t779*t789 &
2200  *phirhob/0.2e1_dp + t795*t789*k_srhob + t801 - t96*t274*( &
2201  -t767*k_frhoa*t523*k_frhob/0.2e1_dp + t266*k_frhoarhob* &
2202  t7)/0.2e1_dp + t1228 + t1231 + t812
2203  arhoarhob = 0.2e1_dp*t820*t822*t432 - t101*t281*( &
2204  -epsilon_c_unifrhoarhob*t100*t105 + 0.3e1_dp*t282*t429 + &
2205  0.3e1_dp*t427*t286 - 0.12e2_dp*t102*t840*phirhob + &
2206  0.3e1_dp*t102*t285*phirhoarhob)*t107 - t851*t289*t432* &
2207  t107
2208  t1269 = t445*t119
2209  t1283 = arhoarhob*t111
2210  t1285 = 0.2e1_dp*t909*trhob
2211  t1286 = arhob*t
2212  t1288 = 0.2e1_dp*t1286*trhoa
2213  t1291 = 0.2e1_dp*a*trhob*trhoa
2214  t1293 = 0.2e1_dp*t303*trhoarhob
2215  t1304 = t445*t312
2216  t1327 = t1283 + t1285 + t1288 + t1291 + t1293 + 0.2e1_dp*arhob* &
2217  t116*arhoa + 0.8e1_dp*t944*arhoa*trhob + 0.2e1_dp*t314* &
2218  arhoarhob + 0.8e1_dp*t944*trhoa*arhob + 0.12e2_dp*t953* &
2219  trhoa*trhob + 0.4e1_dp*t318*trhoarhob
2220  t1330 = 0.2e1_dp*t101*trhob*t299 + 0.2e1_dp*t297*t1269* &
2221  trhoa - 0.2e1_dp*t297*t313*trhoa*t453 + 0.2e1_dp*t297* &
2222  t298*trhoarhob + 0.2e1_dp*t297*t904*trhob + t101*t111*( &
2223  t1283 + t1285 + t1288 + t1291 + t1293)*t119 - t310*t924*t453 - &
2224  0.2e1_dp*t297*t313*t321*trhob - t310*t1304*t321 + &
2225  0.2e1_dp*t310*t936*t321*t453 - t310*t313*t1327
2226 
2227  e_ra_rb(ii) = e_ra_rb(ii) + &
2228  scale_ec*(epsilon_cggarhob + epsilon_cggarhoa + &
2229  my_rho*(epsilon_c_unifrhoarhob + 0.6e1_dp*t858*t294*phirhob + &
2230  0.3e1_dp*t293*t457*phirhoa + 0.3e1_dp*t293*t123* &
2231  phirhoarhob + 0.3e1_dp*t293*t326*phirhob + t110*t1330*t325 &
2232  - t110*t967*t456))
2233 
2234  chirhobrhob = 0.2e1_dp*t163 + 0.2e1_dp*t519
2235  rsrhobrhob = rsrhoarhob
2236  t1342 = t368**2
2237  t1346 = rsrhob**2
2238  e_c_u_0rhobrhob = -0.2e1_dp*t171*rsrhobrhob*t28 + 0.2e1_dp* &
2239  t1033*t1031 - 0.2e1_dp*t543*t1342*t191 + t177*(-t549* &
2240  t1346/0.4e1_dp + t179*rsrhobrhob/0.2e1_dp + beta_2_1* &
2241  rsrhobrhob + 0.3e1_dp/0.4e1_dp*t556*t1346 + 0.3e1_dp/ &
2242  0.2e1_dp*t183*rsrhobrhob + t22*t561*t1346*t564 + t22*t20 &
2243  *rsrhobrhob*t187 - t22*t20*t1346*t564)*t191 + t578* &
2244  t1342*t580*t14/0.2e1_dp
2245  e_c_u_01rhob = e_c_u_0rhob
2246  t1377 = t382**2
2247  t1411 = t396**2
2248  alpha_c1rhob = alpha_crhob
2249  t1440 = chirhob**2
2250  frhobrhob = (0.4e1_dp/0.9e1_dp*t681*t1440 + 0.4e1_dp/ &
2251  0.3e1_dp*t71*chirhobrhob + 0.4e1_dp/0.9e1_dp*t687*t1440 - &
2252  0.4e1_dp/0.3e1_dp*t74*chirhobrhob)*t69
2253  f1rhob = frhob
2254  t1462 = alpha_c1rhob*f
2255  t1465 = alpha_c*f1rhob
2256  t1482 = e_c_u_1rhob - e_c_u_01rhob
2257  t1489 = t1482*f
2258  t1492 = t84*f1rhob
2259  t1501 = -0.4e1_dp*t77*t245*chirhobrhob + (-0.2e1_dp*t194* &
2260  rsrhobrhob*t46 + 0.2e1_dp*t1071*t1069 - 0.2e1_dp*t595* &
2261  t1377*t212 + t200*(-t600*t1346/0.4e1_dp + t201*rsrhobrhob &
2262  /0.2e1_dp + beta_2_2*rsrhobrhob + 0.3e1_dp/0.4e1_dp*t606* &
2263  t1346 + 0.3e1_dp/0.2e1_dp*t205*rsrhobrhob + t40*t611*t1346 &
2264  *t564 + t40*t38*rsrhobrhob*t187 - t40*t38*t1346*t564) &
2265  *t212 + t626*t1377*t628*t34/0.2e1_dp - e_c_u_0rhobrhob)*f &
2266  *t80 + t410*f1rhob*t80 + 0.4e1_dp*t411*t415 + t1482* &
2267  frhob*t80 + t84*frhobrhob*t80 + 0.4e1_dp*t413*t415 + &
2268  0.4e1_dp*t1489*t415 + 0.4e1_dp*t1492*t415 + 0.12e2_dp*t85 &
2269  *t79*t1440 + 0.4e1_dp*t85*t244*chirhobrhob
2270  epsilon_c_unifrhobrhob = e_c_u_0rhobrhob + (0.2e1_dp*t215* &
2271  rsrhobrhob*t64 - 0.2e1_dp*t1106*t1104 + 0.2e1_dp*t643* &
2272  t1411*t233 - t221*(-t648*t1346/0.4e1_dp + t222*rsrhobrhob &
2273  /0.2e1_dp + beta_2_3*rsrhobrhob + 0.3e1_dp/0.4e1_dp*t654* &
2274  t1346 + 0.3e1_dp/0.2e1_dp*t226*rsrhobrhob + t58*t659*t1346 &
2275  *t564 + t58*t56*rsrhobrhob*t187 - t58*t56*t1346*t564) &
2276  *t233 - t674*t1411*t676*t52/0.2e1_dp)*f*t82 + &
2277  alpha_crhob*f1rhob*t82 - 0.4e1_dp*t403*t407 + alpha_c1rhob* &
2278  frhob*t82 + alpha_c*frhobrhob*t82 - 0.4e1_dp*t405*t407 - &
2279  0.4e1_dp*t1462*t407 - 0.4e1_dp*t1465*t407 - 0.12e2_dp*t77 &
2280  *t711*t1440 + t1501
2281  epsilon_c_unif1rhob = e_c_u_01rhob + t1462*t82 + t1465*t82 - &
2282  t409 + t1489*t80 + t1492*t80 + t417
2283  phirhobrhob = -t750*t1440/0.9e1_dp + t257*chirhobrhob/ &
2284  0.3e1_dp - t755*t1440/0.9e1_dp - t259*chirhobrhob/0.3e1_dp
2285  phi1rhob = phirhob
2286  t1514 = k_frhob**2
2287  k_s1rhob = k_srhob
2288  t1520 = t2*phirhob
2289  t1529 = t2*k_srhob
2290  trhobrhob = t775*t1520*phi1rhob + t779*t1520*k_s1rhob/ &
2291  0.2e1_dp + t1228 - t269*t98*phirhobrhob/0.2e1_dp + t779* &
2292  t1529*phi1rhob/0.2e1_dp + t795*t1529*k_s1rhob + t1231 - t96 &
2293  *t274*(-t767*t1514*t523/0.2e1_dp + t266*k_frhoarhob*t7) &
2294  /0.2e1_dp + t269*t277*phi1rhob/0.2e1_dp + t96*t798* &
2295  k_s1rhob/0.2e1_dp + t812
2296  t1rhob = -t269*t98*phi1rhob/0.2e1_dp - t96*t274*k_s1rhob &
2297  /0.2e1_dp - t278/0.2e1_dp
2298  t1550 = epsilon_c_unif1rhob*t100
2299  t1552 = t285*phi1rhob
2300  t1555 = -t1550*t105 + 0.3e1_dp*t102*t1552
2301  arhobrhob = 0.2e1_dp*t820*t432*t821*t1555 - t101*t281* &
2302  (-epsilon_c_unifrhobrhob*t100*t105 + 0.3e1_dp*t427*t1552 + &
2303  0.3e1_dp*t1550*t429 - 0.12e2_dp*t102*t839*phirhob* &
2304  phi1rhob + 0.3e1_dp*t102*t285*phirhobrhob)*t107 - t851* &
2305  t432*t1555*t107
2306  a1rhob = -t101*t281*t1555*t107
2307  t1588 = a1rhob*t111
2308  t1590 = 0.2e1_dp*t303*t1rhob
2309  t1591 = t1588 + t1590
2310  t1599 = t1588 + t1590 + 0.2e1_dp*t314*a1rhob + 0.4e1_dp*t318 &
2311  *t1rhob
2312  t1602 = 0.2e1_dp*t297*t298*t1rhob + t101*t111*t1591* &
2313  t119 - t310*t313*t1599
2314  t1603 = t1602*t325
2315  t1630 = arhobrhob*t111
2316  t1632 = 0.2e1_dp*t1286*t1rhob
2317  t1635 = 0.2e1_dp*a1rhob*t*trhob
2318  t1638 = 0.2e1_dp*a*t1rhob*trhob
2319  t1640 = 0.2e1_dp*t303*trhobrhob
2320  t1674 = t1630 + t1632 + t1635 + t1638 + t1640 + 0.2e1_dp*a1rhob &
2321  *t116*arhob + 0.8e1_dp*t944*arhob*t1rhob + 0.2e1_dp*t314 &
2322  *arhobrhob + 0.8e1_dp*t944*trhob*a1rhob + 0.12e2_dp*t953* &
2323  trhob*t1rhob + 0.4e1_dp*t318*trhobrhob
2324  t1677 = 0.2e1_dp*t101*t1rhob*t439 + 0.2e1_dp*t297*t1591 &
2325  *t119*trhob - 0.2e1_dp*t297*t313*trhob*t1599 + 0.2e1_dp* &
2326  t297*t298*trhobrhob + 0.2e1_dp*t297*t1269*t1rhob + t101* &
2327  t111*(t1630 + t1632 + t1635 + t1638 + t1640)*t119 - t310* &
2328  t1304*t1599 - 0.2e1_dp*t297*t313*t453*t1rhob - t310* &
2329  t1591*t312*t453 + 0.2e1_dp*t310*t936*t453*t1599 - t310* &
2330  t313*t1674
2331  t1680 = t456*t966
2332  kf_brhobrhob = -0.2e1_dp/0.9e1_dp*t124/t460/t142*t763
2333  ex_unif_b1rhob = ex_unif_brhob
2334  t1698 = kf_brhob**2
2335  s_b1rhob = s_brhob
2336  t1711 = 0.1e1_dp/t475/t153*t998
2337  t1712 = t1711*t150
2338  t1713 = s_brhob*t135
2339  fx_b1rhob = 0.2e1_dp*t477*s_b*s_b1rhob
2340 
2341  e_rb_rb(ii) = e_rb_rb(ii) + &
2342  scale_ex*(0.2e1_dp*ex_unif_b1rhob*fx_b + &
2343  0.2e1_dp*ex_unif_b*fx_b1rhob + 0.2e1_dp*ex_unif_brhob*fx_b - &
2344  0.3e1_dp/0.2e1_dp*my_rhob*t7*kf_brhobrhob*fx_b + 0.2e1_dp* &
2345  t481*fx_b1rhob + 0.2e1_dp*ex_unif_b*fx_brhob + 0.2e1_dp*my_rhob &
2346  *ex_unif_b1rhob*fx_brhob + 0.2e1_dp*t156*(-0.8e1_dp*t1712 &
2347  *t1713*s_b1rhob + 0.2e1_dp*t477*s_b1rhob*s_brhob + 0.2e1_dp &
2348  *t477*s_b*(my_norm_drhob/t466/kf_b*t148*t1698 + t468* &
2349  t472*kf_brhob - t468*t148*kf_brhobrhob/0.2e1_dp + t147/ &
2350  t471/my_rhob)))/0.2e1_dp + scale_ec*(epsilon_c_unif1rhob + &
2351  0.3e1_dp*t293*t123*phi1rhob + t110*t1603 + epsilon_cggarhob &
2352  + my_rho*(epsilon_c_unifrhobrhob + 0.6e1_dp*t858*t436*phi1rhob &
2353  + 0.3e1_dp*t293*t1603*phirhob + 0.3e1_dp*t293*t123* &
2354  phirhobrhob + 0.3e1_dp*t293*t457*phi1rhob + t110*t1677* &
2355  t325 - t110*t1680*t1602))
2356  t1739 = t268*t97
2357  t1741 = t95*t273
2358  t1743 = t488*t163
2359  trhoanorm_drho = -t1739*t776/0.2e1_dp - t1741*t789/ &
2360  0.2e1_dp - t1743/0.2e1_dp
2361  t1748 = t101*tnorm_drho
2362  t1765 = t909*tnorm_drho
2363  t1766 = t494*trhoa
2364  t1767 = t303*trhoanorm_drho
2365  t1801 = 0.2e1_dp*t1748*t299 + 0.4e1_dp*t310*t494*t119* &
2366  trhoa - 0.2e1_dp*t297*t313*trhoa*t502 + 0.2e1_dp*t297* &
2367  t298*trhoanorm_drho + 0.2e1_dp*t297*t904*tnorm_drho + t101* &
2368  t111*(0.2e1_dp*t1765 + 0.2e1_dp*t1766 + 0.2e1_dp*t1767)* &
2369  t119 - t310*t924*t502 - 0.2e1_dp*t297*t313*t321* &
2370  tnorm_drho - 0.2e1_dp*t493*t494*t312*t321 + 0.2e1_dp*t310 &
2371  *t936*t321*t502 - t310*t313*(0.2e1_dp*t1765 + 0.2e1_dp* &
2372  t1766 + 0.2e1_dp*t1767 + 0.8e1_dp*t944*arhoa*tnorm_drho + &
2373  0.12e2_dp*t953*trhoa*tnorm_drho + 0.4e1_dp*t318* &
2374  trhoanorm_drho)
2375 
2376  e_ra_ndr(ii) = e_ra_ndr(ii) + &
2377  scale_ec*(hnorm_drho + my_rho*(0.3e1_dp* &
2378  t293*t506*phirhoa + t110*t1801*t325 - t110*t967*t505))
2379 
2380  trhobnorm_drho = -t1739*t1520/0.2e1_dp - t1741*t1529/ &
2381  0.2e1_dp - t1743/0.2e1_dp
2382  t1829 = t1286*tnorm_drho
2383  t1830 = t494*trhob
2384  t1831 = t303*trhobnorm_drho
2385  t1865 = 0.2e1_dp*t1748*t439 + 0.4e1_dp*t310*t494*t119* &
2386  trhob - 0.2e1_dp*t297*t313*trhob*t502 + 0.2e1_dp*t297* &
2387  t298*trhobnorm_drho + 0.2e1_dp*t297*t1269*tnorm_drho + t101 &
2388  *t111*(0.2e1_dp*t1829 + 0.2e1_dp*t1830 + 0.2e1_dp*t1831)* &
2389  t119 - t310*t1304*t502 - 0.2e1_dp*t297*t313*t453* &
2390  tnorm_drho - 0.2e1_dp*t493*t494*t312*t453 + 0.2e1_dp*t310 &
2391  *t936*t453*t502 - t310*t313*(0.2e1_dp*t1829 + 0.2e1_dp* &
2392  t1830 + 0.2e1_dp*t1831 + 0.8e1_dp*t944*arhob*tnorm_drho + &
2393  0.12e2_dp*t953*trhob*tnorm_drho + 0.4e1_dp*t318* &
2394  trhobnorm_drho)
2395 
2396  e_rb_ndr(ii) = e_rb_ndr(ii) + &
2397  scale_ec*(hnorm_drho + my_rho*(0.3e1_dp* &
2398  t293*t506*phirhob + t110*t1865*t325 - t110*t1680*t505))
2399 
2400  t1871 = tnorm_drho**2
2401  t1876 = a*t1871
2402  t1888 = t502**2
2403  t1901 = t505**2
2404 
2405  e_ndr_ndr(ii) = e_ndr_ndr(ii) + &
2406  scale_ec*my_rho*(t110*(0.2e1_dp* &
2407  t101*t1871*t113*t119 + 0.10e2_dp*t310*t1876*t119 - &
2408  0.4e1_dp*t297*t313*tnorm_drho*t502 - 0.4e1_dp*t493*t494 &
2409  *t312*t502 + 0.2e1_dp*t310*t936*t1888 - t310*t313*( &
2410  0.2e1_dp*t1876 + 0.12e2_dp*t953*t1871))*t325 - t110*t1901 &
2411  *t966)
2412 
2413  e_ra_ndra(ii) = e_ra_ndra(ii) + &
2414  scale_ex*(0.2e1_dp*ex_unif_a* &
2415  fx_anorm_drhoa + 0.2e1_dp*t350*fx_anorm_drhoa + 0.2e1_dp*t140 &
2416  *(-0.8e1_dp*t1000*t1001*s_anorm_drhoa + 0.2e1_dp*t346* &
2417  s_anorm_drhoa*s_arhoa + 0.2e1_dp*t346*s_a*(-t336*t131* &
2418  kf_arhoa/0.2e1_dp - t129*t341/0.2e1_dp)))/0.2e1_dp
2419 
2420  t1922 = s_anorm_drhoa**2
2421 
2422  e_ndra_ndra(ii) = e_ndra_ndra(ii) + &
2423  scale_ex*t140*(-0.8e1_dp*t999* &
2424  t133*t1922*t135 + 0.2e1_dp*t346*t1922)
2425 
2426  e_rb_ndrb(ii) = e_rb_ndrb(ii) + &
2427  scale_ex*(0.2e1_dp*ex_unif_b* &
2428  fx_bnorm_drhob + 0.2e1_dp*t481*fx_bnorm_drhob + 0.2e1_dp*t156 &
2429  *(-0.8e1_dp*t1712*t1713*s_bnorm_drhob + 0.2e1_dp*t477* &
2430  s_bnorm_drhob*s_brhob + 0.2e1_dp*t477*s_b*(-t467*t148* &
2431  kf_brhob/0.2e1_dp - t146*t472/0.2e1_dp)))/0.2e1_dp
2432 
2433  t1949 = s_bnorm_drhob**2
2434  e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii) + &
2435  scale_ex*t156*(-0.8e1_dp*t1711* &
2436  t150*t1949*t135 + 0.2e1_dp*t477*t1949)
2437  END IF
2438  END IF
2439  END DO
2440 !$OMP END DO
2441  END SELECT
2442  END SUBROUTINE pbe_lsd_calc
2443 
2444 END MODULE xc_pbe
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::perdew1996
integer, save, public perdew1996
Definition: bibliography.F:43
bibliography::perdew2008
integer, save, public perdew2008
Definition: bibliography.F:43
bibliography::zhang1998
integer, save, public zhang1998
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
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
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_norm_drho
integer, parameter, public deriv_norm_drho
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition: xc_derivative_desc.F:35
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_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
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_input_constants
input constants for xc
Definition: xc_input_constants.F:10
xc_input_constants::xc_pbe_orig
integer, parameter, public xc_pbe_orig
Definition: xc_input_constants.F:33
xc_input_constants::xc_pbe_rev
integer, parameter, public xc_pbe_rev
Definition: xc_input_constants.F:33
xc_input_constants::xc_pbe_sol
integer, parameter, public xc_pbe_sol
Definition: xc_input_constants.F:33
xc_pbe
calculates the pbe correlation functional
Definition: xc_pbe.F:21
xc_pbe::pbe_lda_eval
subroutine, public pbe_lda_eval(rho_set, deriv_set, grad_deriv, pbe_params)
evaluates the pbe correlation functional for lda
Definition: xc_pbe.F:280
xc_pbe::pbe_lsd_info
subroutine, public pbe_lsd_info(pbe_params, reference, shortform, needs, max_deriv)
return various information on the functional
Definition: xc_pbe.F:174
xc_pbe::pbe_lsd_eval
subroutine, public pbe_lsd_eval(rho_set, deriv_set, grad_deriv, pbe_params)
evaluates the becke 88 exchange functional for lsd
Definition: xc_pbe.F:1326
xc_pbe::pbe_lda_info
subroutine, public pbe_lda_info(pbe_params, reference, shortform, needs, max_deriv)
return various information on the functional
Definition: xc_pbe.F:70
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