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