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