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