(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
21MODULE 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
57CONTAINS
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
2444END 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
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:49
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78