(git:c641773)
Loading...
Searching...
No Matches
xc_b97.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief calculates the b97 correlation functional
10!> \note
11!> This was generated with the help of a maple worksheet that you can
12!> find under doc/b97.mw .
13!> I did not add 3. derivatives in the polarized (lsd) case because it
14!> would have added lots 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!> Other parametrizations (B97-1 by FA Hamprecht, AJ Cohen, DJ Tozer, NC Handy)
18!> could be easily added, the maple file should be straightforward to extend
19!> to HCTH (and thus implement it also for unrestricted calculations).
20!> \par History
21!> 01.2009 created [fawzi]
22!> \author fawzi
23! **************************************************************************************************
24MODULE xc_b97
25 USE bibliography, ONLY: becke1997,&
26 grimme2006,&
27 cite_reference
28 USE input_section_types, ONLY: section_vals_type,&
29 section_vals_val_get
30 USE kinds, ONLY: dp
31 USE mathconstants, ONLY: pi
32 USE xc_derivative_desc, ONLY: deriv_norm_drho,&
33 deriv_norm_drhoa,&
34 deriv_norm_drhob,&
35 deriv_rho,&
36 deriv_rhoa,&
37 deriv_rhob
38 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
39 xc_dset_get_derivative
40 USE xc_derivative_types, ONLY: xc_derivative_get,&
41 xc_derivative_type
42 USE xc_input_constants, ONLY: xc_b97_3c,&
43 xc_b97_grimme,&
44 xc_b97_mardirossian,&
45 xc_b97_orig
46 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
47 USE xc_rho_set_types, ONLY: xc_rho_set_get,&
48 xc_rho_set_type
49#include "../base/base_uses.f90"
50
51 IMPLICIT NONE
52 PRIVATE
53
54 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_b97'
56
57 PUBLIC :: b97_lda_info, b97_lsd_info, b97_lda_eval, b97_lsd_eval
58
59 REAL(dp), DIMENSION(10) :: params_b97_orig = [0.8094_dp, 0.5073_dp, 0.7481_dp, &
60 0.9454_dp, 0.7471_dp, -4.5961_dp, 0.1737_dp, 2.3487_dp, -2.4868_dp, 1.0_dp - 0.1943_dp]
61 REAL(dp), DIMENSION(10) :: params_b97_grimme = [1.08662_dp, -0.52127_dp, 3.25429_dp, &
62 0.69041_dp, 6.30270_dp, -14.9712_dp, 0.22340_dp, -1.56208_dp, 1.94293_dp, 1.0_dp]
63 REAL(dp), DIMENSION(10) :: params_b97_mardirossian = [0.833_dp, 0.603_dp, 1.194_dp, &
64 1.219_dp, -1.850_dp, 0.00_dp, 0.556_dp, -0.257_dp, 0.00_dp, 1.0_dp]
65 REAL(dp), DIMENSION(10) :: params_b97_3c = [1.076616_dp, -0.469912_dp, 3.322442_dp, &
66 0.635047_dp, 5.532103_dp, -15.301575_dp, 0.543788_dp, -1.444420_dp, 1.637436_dp, 1.0_dp]
67
68CONTAINS
69
70! **************************************************************************************************
71!> \brief ...
72!> \param param ...
73!> \param lda ...
74!> \param sx ...
75!> \param sc ...
76!> \param reference ...
77!> \param shortform ...
78! **************************************************************************************************
79 SUBROUTINE b97_ref(param, lda, sx, sc, reference, shortform)
80 INTEGER :: param
81 LOGICAL :: lda
82 REAL(dp), INTENT(in) :: sx, sc
83 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
84
85 CHARACTER(20) :: pol
86
87 IF (lda) THEN
88 pol = "(unpolarized)"
89 ELSE
90 pol = "(polarized)"
91 END IF
92 SELECT CASE (param)
93 CASE (xc_b97_orig)
94 CALL cite_reference(becke1997)
95 IF (sx == 1._dp .AND. sc == 1._dp) THEN
96 IF (PRESENT(reference)) THEN
97 reference = "A.D.Becke, "// &
98 "J. Chem. Phys, vol. 107, pp. 8554, (1997), needs exact x, "// &
99 pol
100 END IF
101 IF (PRESENT(shortform)) THEN
102 shortform = "B97, Becke 1997 xc functional, needs exact x "//pol
103 END IF
104 ELSE
105 IF (PRESENT(reference)) THEN
106 WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a,a)") &
107 "A.D.Becke, ", &
108 "J. Chem. Phys, vol. 107, pp. 8554, (1997)", &
109 sx, sc, ", needs exact x ", pol
110 END IF
111 IF (PRESENT(shortform)) THEN
112 WRITE (shortform, "(a,a,'sx=',f5.3,'sc=',f5.3)") &
113 "B97, Becke 1997 xc functional (unpolarized)", pol, sx, sc
114 END IF
115 END IF
116 CASE (xc_b97_grimme)
117 CALL cite_reference(becke1997)
118 CALL cite_reference(grimme2006)
119 IF (sx == 1._dp .AND. sc == 1._dp) THEN
120 IF (PRESENT(reference)) THEN
121 reference = "B97-D, Grimme xc functional,"// &
122 " J Comput Chem 27: 1787-1799 (2006),"// &
123 " needs C6 dispersion "//pol
124 END IF
125 IF (PRESENT(shortform)) THEN
126 shortform = "B97-D, Grimme xc functional "//pol
127 END IF
128 ELSE
129 IF (PRESENT(reference)) THEN
130 WRITE (reference, "(a,a,3x,' sx=',f6.3,' sc=',f6.3,1x,a)") &
131 "B97-D, Grimme xc functional,", &
132 " J Comput Chem 27: 1787-1799 (2006) ", &
133 sx, sc, pol
134 END IF
135 IF (PRESENT(shortform)) THEN
136 WRITE (shortform, "(a,a,3x,' sx=',f6.3,' sc=',f6.3,' (LDA)')") &
137 "B97-D, Grimme xc functional ", pol, sx, sc
138 END IF
139 END IF
140 CASE (xc_b97_mardirossian)
141 CALL cite_reference(becke1997)
142! CALL cite_reference(Mardirossian2014)
143 IF (sx == 1._dp .AND. sc == 1._dp) THEN
144 IF (PRESENT(reference)) THEN
145 reference = "wB97X-V, xc functional,"// &
146 " Mardirossian and Head-Gordon, PCCP DOI: 10.1039/c3cp54374a,"// &
147 " needs HFX exchange and VV10 dispersion (NOT TESTED)"//pol
148 END IF
149 IF (PRESENT(shortform)) THEN
150 shortform = "wB97X-V, HFX+B97+VV10 functional (NOT TESTED)"//pol
151 END IF
152 ELSE
153 cpabort("Unsupported scaling factors")
154 END IF
155 CASE (xc_b97_3c)
156 CALL cite_reference(becke1997)
157! CALL cite_reference(Mardirossian2014)
158 IF (sx == 1._dp .AND. sc == 1._dp) THEN
159 IF (PRESENT(reference)) THEN
160 reference = "B97-3c composite method,"// &
161 " S. Grimme, A. Hansen, no reference available, yet,"// &
162 " use TZVP basis set, D3(BJ) dispersion correction"// &
163 " with CALCULATE_C9_TERM and SRB correction"//pol
164 END IF
165 IF (PRESENT(shortform)) THEN
166 shortform = "B97-3c composite method"//pol
167 END IF
168 ELSE
169 cpabort("Unsupported scaling factors")
170 END IF
171 CASE default
172 cpabort("Unsupported parametrization")
173 END SELECT
174 END SUBROUTINE b97_ref
175
176! **************************************************************************************************
177!> \brief return various information on the functional
178!> \param b97_params section selecting the various parameters for the functional
179!> \param reference string with the reference of the actual functional
180!> \param shortform string with the shortform of the functional name
181!> \param needs the components needed by this functional are set to
182!> true (does not set the unneeded components to false)
183!> \param max_deriv ...
184!> \author fawzi
185! **************************************************************************************************
186 SUBROUTINE b97_lda_info(b97_params, reference, shortform, needs, max_deriv)
187 TYPE(section_vals_type), POINTER :: b97_params
188 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
189 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
190 INTEGER, INTENT(out), OPTIONAL :: max_deriv
191
192 INTEGER :: param
193 REAL(kind=dp) :: sc, sx
194
195 CALL section_vals_val_get(b97_params, "parametrization", i_val=param)
196 CALL section_vals_val_get(b97_params, "scale_x", r_val=sx)
197 CALL section_vals_val_get(b97_params, "scale_c", r_val=sc)
198
199 CALL b97_ref(param, .true., sx, sc, reference, shortform)
200 IF (PRESENT(needs)) THEN
201 needs%rho = .true.
202 needs%norm_drho = .true.
203 END IF
204 IF (PRESENT(max_deriv)) max_deriv = 2
205
206 END SUBROUTINE b97_lda_info
207
208! **************************************************************************************************
209!> \brief return various information on the functional
210!> \param b97_params section selecting the various parameters for the functional
211!> \param reference string with the reference of the actual functional
212!> \param shortform string with the shortform of the functional name
213!> \param needs the components needed by this functional are set to
214!> true (does not set the unneeded components to false)
215!> \param max_deriv ...
216!> \author fawzi
217! **************************************************************************************************
218 SUBROUTINE b97_lsd_info(b97_params, reference, shortform, needs, max_deriv)
219 TYPE(section_vals_type), POINTER :: b97_params
220 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
221 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
222 INTEGER, INTENT(out), OPTIONAL :: max_deriv
223
224 INTEGER :: param
225 REAL(kind=dp) :: sc, sx
226
227 CALL section_vals_val_get(b97_params, "parametrization", i_val=param)
228 CALL section_vals_val_get(b97_params, "scale_x", r_val=sx)
229 CALL section_vals_val_get(b97_params, "scale_c", r_val=sc)
230
231 CALL b97_ref(param, .false., sx, sc, reference, shortform)
232 IF (PRESENT(needs)) THEN
233 needs%rho_spin = .true.
234 needs%norm_drho_spin = .true.
235 END IF
236 IF (PRESENT(max_deriv)) max_deriv = 2
237
238 END SUBROUTINE b97_lsd_info
239
240! **************************************************************************************************
241!> \brief evaluates the b97 correlation functional for lda
242!> \param rho_set the density where you want to evaluate the functional
243!> \param deriv_set place where to store the functional derivatives (they are
244!> added to the derivatives)
245!> \param grad_deriv maximum degree of the derivative that should be evaluated
246!> \param b97_params ...
247!> \author fawzi
248! **************************************************************************************************
249 SUBROUTINE b97_lda_eval(rho_set, deriv_set, grad_deriv, b97_params)
250 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
251 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
252 INTEGER, INTENT(in) :: grad_deriv
253 TYPE(section_vals_type), POINTER :: b97_params
254
255 CHARACTER(len=*), PARAMETER :: routinen = 'b97_lda_eval'
256
257 INTEGER :: handle, npoints, param
258 INTEGER, DIMENSION(2, 3) :: bo
259 REAL(kind=dp) :: epsilon_norm_drho, epsilon_rho, scale_c, &
260 scale_x
261 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
262 e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
263 e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
264 TYPE(xc_derivative_type), POINTER :: deriv
265
266 CALL timeset(routinen, handle)
267
268 CALL xc_rho_set_get(rho_set, rho=rho, &
269 norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
270 drho_cutoff=epsilon_norm_drho)
271 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
272
273 dummy => rho
274
275 e_0 => dummy
276 e_rho => dummy
277 e_ndrho => dummy
278 e_rho_rho => dummy
279 e_ndrho_rho => dummy
280 e_ndrho_ndrho => dummy
281 e_rho_rho_rho => dummy
282 e_ndrho_rho_rho => dummy
283 e_ndrho_ndrho_rho => dummy
284 e_ndrho_ndrho_ndrho => dummy
285
286 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
287 allocate_deriv=.true.)
288 CALL xc_derivative_get(deriv, deriv_data=e_0)
289 IF (grad_deriv >= 1) THEN
290 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
291 allocate_deriv=.true.)
292 CALL xc_derivative_get(deriv, deriv_data=e_rho)
293 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
294 allocate_deriv=.true.)
295 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
296 END IF
297 IF (grad_deriv >= 2) THEN
298 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
299 allocate_deriv=.true.)
300 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
301 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
302 allocate_deriv=.true.)
303 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
304 deriv => xc_dset_get_derivative(deriv_set, &
305 [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
306 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
307 END IF
308 IF (grad_deriv >= 3) THEN
309 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
310 allocate_deriv=.true.)
311 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
312 deriv => xc_dset_get_derivative(deriv_set, &
313 [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.true.)
314 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
315 deriv => xc_dset_get_derivative(deriv_set, &
316 [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.true.)
317 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
318 deriv => xc_dset_get_derivative(deriv_set, &
319 [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
320 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
321 END IF
322 IF (grad_deriv > 3) THEN
323 cpabort("derivatives bigger than 3 not implemented")
324 END IF
325
326 CALL section_vals_val_get(b97_params, "parametrization", i_val=param)
327 CALL section_vals_val_get(b97_params, "scale_c", r_val=scale_c)
328 CALL section_vals_val_get(b97_params, "scale_x", r_val=scale_x)
329
330!$OMP PARALLEL DEFAULT(NONE) &
331!$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho, e_rho_rho) &
332!$OMP SHARED(e_ndrho_rho, e_ndrho_ndrho) &
333!$OMP SHARED(grad_deriv, npoints, epsilon_rho) &
334!$OMP SHARED(epsilon_norm_drho, param, scale_c, scale_x)
335
336 CALL b97_lda_calc(rho_tot=rho, norm_drho=norm_drho, &
337 e_0=e_0, e_r=e_rho, e_ndr=e_ndrho, e_r_r=e_rho_rho, &
338 e_r_ndr=e_ndrho_rho, e_ndr_ndr=e_ndrho_ndrho, &
339 ! e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho,&
340 ! e_ndrho_ndrho_rho=e_ndrho_ndrho_rho,e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho,&
341 grad_deriv=grad_deriv, &
342 npoints=npoints, epsilon_rho=epsilon_rho, &
343 param=param, scale_c_in=scale_c, scale_x_in=scale_x)
344
345!$OMP END PARALLEL
346
347 NULLIFY (dummy)
348 CALL timestop(handle)
349 END SUBROUTINE b97_lda_eval
350
351! **************************************************************************************************
352!> \brief evaluates the b 97 xc functional for lsd
353!> \param rho_set ...
354!> \param deriv_set ...
355!> \param grad_deriv ...
356!> \param b97_params ...
357!> \author fawzi
358! **************************************************************************************************
359 SUBROUTINE b97_lsd_eval(rho_set, deriv_set, grad_deriv, b97_params)
360 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
361 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
362 INTEGER, INTENT(in) :: grad_deriv
363 TYPE(section_vals_type), POINTER :: b97_params
364
365 CHARACTER(len=*), PARAMETER :: routinen = 'b97_lsd_eval'
366
367 INTEGER :: handle, npoints, param
368 INTEGER, DIMENSION(2, 3) :: bo
369 REAL(kind=dp) :: epsilon_rho, scale_c, scale_x
370 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndra, e_ndra_ndra, &
371 e_ndra_ndrb, e_ndra_ra, e_ndra_rb, e_ndrb, e_ndrb_ndrb, e_ndrb_ra, e_ndrb_rb, e_ra, &
372 e_ra_ra, e_ra_rb, e_rb, e_rb_rb, norm_drhoa, norm_drhob, rhoa, rhob
373 TYPE(xc_derivative_type), POINTER :: deriv
374
375 CALL timeset(routinen, handle)
376 NULLIFY (deriv)
377
378 CALL xc_rho_set_get(rho_set, &
379 rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
380 norm_drhob=norm_drhob, &
381 rho_cutoff=epsilon_rho, &
382 local_bounds=bo)
383 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
384
385 dummy => rhoa
386 e_0 => dummy
387 e_ra => dummy
388 e_rb => dummy
389 e_ndra => dummy
390 e_ndrb => dummy
391
392 e_ndra_ra => dummy
393 e_ndra_rb => dummy
394 e_ndrb_rb => dummy
395 e_ndrb_ra => dummy
396
397 e_ndra_ndra => dummy
398 e_ndrb_ndrb => dummy
399 e_ndra_ndrb => dummy
400
401 e_ra_ra => dummy
402 e_ra_rb => dummy
403 e_rb_rb => dummy
404
405 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
406 allocate_deriv=.true.)
407 CALL xc_derivative_get(deriv, deriv_data=e_0)
408 IF (grad_deriv >= 1) THEN
409 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
410 allocate_deriv=.true.)
411 CALL xc_derivative_get(deriv, deriv_data=e_ra)
412 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
413 allocate_deriv=.true.)
414 CALL xc_derivative_get(deriv, deriv_data=e_rb)
415 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
416 allocate_deriv=.true.)
417 CALL xc_derivative_get(deriv, deriv_data=e_ndra)
418 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
419 allocate_deriv=.true.)
420 CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
421 END IF
422 IF (grad_deriv >= 2) THEN
423 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
424 allocate_deriv=.true.)
425 CALL xc_derivative_get(deriv, deriv_data=e_ra_ra)
426 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
427 allocate_deriv=.true.)
428 CALL xc_derivative_get(deriv, deriv_data=e_ra_rb)
429 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
430 allocate_deriv=.true.)
431 CALL xc_derivative_get(deriv, deriv_data=e_rb_rb)
432 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
433 allocate_deriv=.true.)
434 CALL xc_derivative_get(deriv, deriv_data=e_ndra_ra)
435 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
436 allocate_deriv=.true.)
437 CALL xc_derivative_get(deriv, deriv_data=e_ndra_rb)
438 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
439 allocate_deriv=.true.)
440 CALL xc_derivative_get(deriv, deriv_data=e_ndrb_rb)
441 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa], &
442 allocate_deriv=.true.)
443 CALL xc_derivative_get(deriv, deriv_data=e_ndrb_ra)
444 deriv => xc_dset_get_derivative(deriv_set, &
445 [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.true.)
446 CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndra)
447 deriv => xc_dset_get_derivative(deriv_set, &
448 [deriv_norm_drhoa, deriv_norm_drhob], allocate_deriv=.true.)
449 CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndrb)
450 deriv => xc_dset_get_derivative(deriv_set, &
451 [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.true.)
452 CALL xc_derivative_get(deriv, deriv_data=e_ndrb_ndrb)
453 END IF
454 IF (grad_deriv >= 3) THEN
455 cpabort("Derivatives higher than 3 not implemented!")
456 END IF
457
458 CALL section_vals_val_get(b97_params, "parametrization", i_val=param)
459 CALL section_vals_val_get(b97_params, "scale_x", r_val=scale_x)
460 CALL section_vals_val_get(b97_params, "scale_c", r_val=scale_c)
461
462!$OMP PARALLEL DEFAULT (NONE) &
463!$OMP SHARED(rhoa, rhob, norm_drhoa, norm_drhob, e_0, e_ra) &
464!$OMP SHARED(e_rb, e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb) &
465!$OMP SHARED(e_ndra_ra, e_ndrb_ra, e_ndrb_rb, e_ndra_rb) &
466!$OMP SHARED(e_ndra_ndra, e_ndrb_ndrb, e_ndra_ndrb) &
467!$OMP SHARED(grad_deriv, npoints, param, scale_c, scale_x) &
468!$OMP SHARED(epsilon_rho)
469
470 CALL b97_lsd_calc( &
471 rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
472 norm_drhob=norm_drhob, e_0=e_0, &
473 e_ra=e_ra, e_rb=e_rb, &
474 e_ndra=e_ndra, e_ndrb=e_ndrb, &
475 e_ra_ra=e_ra_ra, e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, &
476 e_ra_ndra=e_ndra_ra, e_ra_ndrb=e_ndrb_ra, &
477 e_rb_ndrb=e_ndrb_rb, e_rb_ndra=e_ndra_rb, &
478 e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, &
479 e_ndra_ndrb=e_ndra_ndrb, &
480 grad_deriv=grad_deriv, npoints=npoints, &
481 epsilon_rho=epsilon_rho, &
482 param=param, scale_c_in=scale_c, scale_x_in=scale_x)
483
484!$OMP END PARALLEL
485
486 NULLIFY (dummy)
487 CALL timestop(handle)
488 END SUBROUTINE b97_lsd_eval
489
490! **************************************************************************************************
491!> \brief ...
492!> \param param ...
493!> \return ...
494! **************************************************************************************************
495 FUNCTION b97_coeffs(param) RESULT(res)
496 INTEGER, INTENT(in) :: param
497 REAL(dp), DIMENSION(10) :: res
498
499 SELECT CASE (param)
500 CASE (xc_b97_orig)
501 res = params_b97_orig
502 CASE (xc_b97_grimme)
503 res = params_b97_grimme
504 CASE (xc_b97_mardirossian)
505 res = params_b97_mardirossian
506 CASE (xc_b97_3c)
507 res = params_b97_3c
508 CASE default
509 cpabort("Unsupported parametrization")
510 res = 0.0_dp
511 END SELECT
512 END FUNCTION b97_coeffs
513
514! **************************************************************************************************
515!> \brief low level calculation of the b97 exchange-correlation functional for lsd
516!> \param rhoa alpha spin density
517!> \param rhob beta spin desnity
518!> \param norm_drhoa || grad rhoa ||
519!> \param norm_drhob || grad rhoa ||
520!> \param e_0 adds to it the local value of the functional
521!> \param e_ra derivative of the functional with respect to ra
522!> \param e_rb derivative of the functional with respect to rb
523!> \param e_ndra derivative of the functional with respect to ndra
524!> \param e_ndrb derivative of the functional with respect to ndrb
525!> \param e_ra_ndra derivative of the functional with respect to ra_ndra
526!> \param e_ra_ndrb derivative of the functional with respect to ra_ndrb
527!> \param e_rb_ndra derivative of the functional with respect to rb_ndra
528!> \param e_rb_ndrb derivative of the functional with respect to rb_ndrb
529!> \param e_ndra_ndra derivative of the functional with respect to ndra_ndra
530!> \param e_ndrb_ndrb derivative of the functional with respect to ndrb_ndrb
531!> \param e_ndra_ndrb derivative of the functional with respect to ndra_ndrb
532!> \param e_ra_ra derivative of the functional with respect to ra_ra
533!> \param e_ra_rb derivative of the functional with respect to ra_rb
534!> \param e_rb_rb derivative of the functional with respect to rb_rb
535!> \param grad_deriv ...
536!> \param npoints ...
537!> \param epsilon_rho ...
538!> \param param ...
539!> \param scale_c_in derivative of the functional with respect to c_in
540!> \param scale_x_in derivative of the functional with respect to x_in
541!> \author fawzi
542! **************************************************************************************************
543 SUBROUTINE b97_lsd_calc(rhoa, rhob, norm_drhoa, norm_drhob, &
544 e_0, e_ra, e_rb, e_ndra, e_ndrb, &
545 e_ra_ndra, e_ra_ndrb, e_rb_ndra, e_rb_ndrb, &
546 e_ndra_ndra, e_ndrb_ndrb, e_ndra_ndrb, &
547 e_ra_ra, e_ra_rb, e_rb_rb, &
548 grad_deriv, npoints, epsilon_rho, &
549 param, scale_c_in, scale_x_in)
550 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drhoa, norm_drhob
551 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ndra, e_ndrb, e_ra_ndra, &
552 e_ra_ndrb, e_rb_ndra, e_rb_ndrb, e_ndra_ndra, e_ndrb_ndrb, e_ndra_ndrb, e_ra_ra, e_ra_rb, &
553 e_rb_rb
554 INTEGER, INTENT(in) :: grad_deriv, npoints
555 REAL(kind=dp), INTENT(in) :: epsilon_rho
556 INTEGER, INTENT(in) :: param
557 REAL(kind=dp), INTENT(in) :: scale_c_in, scale_x_in
558
559 INTEGER :: ii
560 REAL(kind=dp) :: a_1, a_2, a_3, alpha_1_1, alpha_1_2, alpha_1_3, alpha_c, alpha_c1rhoa, &
561 alpha_c1rhob, alpha_crhoa, alpha_crhob, beta_1_1, beta_1_2, beta_1_3, beta_2_1, beta_2_2, &
562 beta_2_3, beta_3_1, beta_3_2, beta_3_3, beta_4_1, beta_4_2, beta_4_3, c_cab_0, c_cab_1, &
563 c_cab_2, c_css_0, c_css_1, c_css_2, c_x_0, c_x_1, c_x_2, chi, chirhoa, chirhoarhoa, &
564 chirhoarhob, chirhob, chirhobrhob, e_c_u_0, e_c_u_01rhoa, e_c_u_01rhob, e_c_u_0rhoa, &
565 e_c_u_0rhoarhoa, e_c_u_0rhoarhob, e_c_u_0rhob, e_c_u_0rhobrhob, e_c_u_1rhoa, e_c_u_1rhob, &
566 e_lsda_c_a, e_lsda_c_a1rhoa, e_lsda_c_ab, e_lsda_c_abrhoa
567 REAL(kind=dp) :: e_lsda_c_abrhob, e_lsda_c_arhoa, e_lsda_c_arhoarhoa, e_lsda_c_b, &
568 e_lsda_c_b1rhob, e_lsda_c_brhob, e_lsda_c_brhobrhob, e_lsda_x_a, e_lsda_x_arhoa, &
569 e_lsda_x_b, e_lsda_x_brhob, epsilon_c_unif, epsilon_c_unif1rhoa, epsilon_c_unif1rhob, &
570 epsilon_c_unif_a, epsilon_c_unif_a1rhoa, epsilon_c_unif_arhoa, epsilon_c_unif_b, &
571 epsilon_c_unif_b1rhob, epsilon_c_unif_brhob, epsilon_c_unifrhoa, epsilon_c_unifrhob, exc, &
572 exc_norm_drhoa, exc_norm_drhoa_norm_drhoa, exc_norm_drhoa_norm_drhob, exc_norm_drhob, &
573 exc_norm_drhob_norm_drhob, exc_rhoa, exc_rhoa_norm_drhoa, exc_rhoa_norm_drhob
574 REAL(kind=dp) :: exc_rhoa_rhoa, exc_rhoa_rhob, exc_rhob, exc_rhob_norm_drhoa, &
575 exc_rhob_norm_drhob, exc_rhob_rhob, f, f1rhoa, f1rhob, frhoa, frhoarhoa, frhoarhob, &
576 frhob, frhobrhob, gamma_c_ab, gamma_c_ss, gamma_x, gc_a, gc_ab, gc_abnorm_drhoa, &
577 gc_abnorm_drhob, gc_abrhoa, gc_abrhob, gc_anorm_drhoa, gc_arhoa, gc_b, gc_bnorm_drhob, &
578 gc_brhob, gx_a, gx_anorm_drhoa, gx_arhoa, gx_b, gx_bnorm_drhob, gx_brhob, my_norm_drhoa, &
579 my_norm_drhob, my_rhoa, my_rhob, p_1, p_2, p_3, rho, rs, rs_a, rs_arhoa, rs_arhoarhoa, &
580 rs_b, rs_brhob, rs_brhobrhob, rsrhoa, rsrhoarhoa, rsrhoarhob, rsrhob, rsrhobrhob, s_a
581 REAL(kind=dp) :: s_a_2, s_a_21norm_drhoa, s_a_21rhoa, s_a_2norm_drhoa, &
582 s_a_2norm_drhoanorm_drhoa, s_a_2rhoa, s_a_2rhoanorm_drhoa, s_a_2rhoarhoa, s_anorm_drhoa, &
583 s_arhoa, s_arhoanorm_drhoa, s_arhoarhoa, s_avg_2, s_avg_21norm_drhoa, s_avg_21norm_drhob, &
584 s_avg_21rhoa, s_avg_21rhob, s_avg_2norm_drhoa, s_avg_2norm_drhoanorm_drhoa, &
585 s_avg_2norm_drhob, s_avg_2norm_drhobnorm_drhob, s_avg_2rhoa, s_avg_2rhoanorm_drhoa, &
586 s_avg_2rhoarhoa, s_avg_2rhob, s_avg_2rhobnorm_drhob, s_avg_2rhobrhob, s_b, s_b_2, &
587 s_b_21norm_drhob, s_b_21rhob, s_b_2norm_drhob, s_b_2norm_drhobnorm_drhob, s_b_2rhob, &
588 s_b_2rhobnorm_drhob
589 REAL(kind=dp) :: s_b_2rhobrhob, s_bnorm_drhob, s_brhob, s_brhobnorm_drhob, s_brhobrhob, &
590 scale_c, scale_x, t1, t101, t1012, t1014, t102, t1047, t1049, t105, t106, t107, t108, &
591 t110, t1107, t112, t113, t1136, t1152, t1157, t116, t1161, t1165, t117, t1172, t1173, &
592 t1175, t12, t120, t1201, t1205, t122, t1236, t125, t127, t1270, t128, t129, t1299, t13, &
593 t1321, t1324, t133, t134, t1341, t1348, t1351, t1360, t1368, t138, t1388, t139, t1394, &
594 t14, t1406, t1407, t1419, t142, t1422, t1436, t1437, t144, t1440, t1451, t1452, t1455, &
595 t1457, t1458, t147, t149, t15, t150, t151, t155, t156, t1571, t1589, t1590
596 REAL(kind=dp) :: t1593, t16, t160, t1605, t161, t162, t163, t164, t165, t166, t167, t168, &
597 t170, t1719, t173, t1738, t1753, t176, t18, t186, t188, t191, t192, t194, t196, t197, &
598 t198, t199, t2, t20, t207, t208, t209, t21, t210, t212, t216, t220, t221, t222, t223, &
599 t224, t228, t232, t235, t236, t237, t239, t243, t244, t245, t246, t25, t250, t256, t257, &
600 t258, t26, t260, t264, t265, t266, t267, t27, t271, t277, t278, t279, t28, t285, t287, &
601 t289, t29, t290, t291, t293, t294, t295, t297, t299, t3, t301, t302, t304, t31, t312, &
602 t313, t314, t315, t316, t320, t324, t327, t328, t33, t336, t337, t338
603 REAL(kind=dp) :: t339, t34, t344, t345, t346, t347, t35, t36, t365, t367, t37, t370, t371, &
604 t374, t375, t376, t377, t396, t4, t40, t410, t42, t424, t43, t431, t433, t435, t437, &
605 t438, t439, t441, t443, t445, t446, t448, t456, t457, t458, t459, t46, t460, t464, t468, &
606 t471, t472, t48, t480, t484, t485, t486, t49, t50, t51, t512, t516, t539, t543, t55, &
607 t555, t56, t560, t564, t568, t575, t576, t577, t579, t6, t60, t600, t601, t605, t606, &
608 t608, t619, t62, t621, t626, t627, t631, t632, t633, t639, t644, t646, t647, t659, t66, &
609 t661, t662, t663, t67, t671, t673, t678, t679, t68, t683, t689, t69
610 REAL(kind=dp) :: t694, t7, t707, t709, t710, t711, t719, t721, t726, t727, t73, t731, t737, &
611 t74, t742, t755, t757, t758, t759, t764, t765, t766, t771, t772, t78, t790, t793, t796, &
612 t8, t80, t811, t818, t821, t830, t838, t84, t85, t858, t86, t864, t87, t876, t877, t889, &
613 t892, t906, t907, t91, t911, t913, t914, t92, t925, t926, t929, t930, t932, t933, t94, &
614 t97, t974, t976, t98, t981, t99, t993, u_c_a, u_c_a1rhoa, u_c_ab, u_c_ab1rhoa, &
615 u_c_ab1rhob, u_c_abnorm_drhoa, u_c_abnorm_drhoanorm_drhoa, u_c_abnorm_drhoanorm_drhob, &
616 u_c_abnorm_drhob, u_c_abnorm_drhobnorm_drhob, u_c_abrhoa
617 REAL(kind=dp) :: u_c_abrhoanorm_drhoa, u_c_abrhoanorm_drhob, u_c_abrhoarhoa, u_c_abrhoarhob, &
618 u_c_abrhob, u_c_abrhobnorm_drhoa, u_c_abrhobnorm_drhob, u_c_abrhobrhob, u_c_anorm_drhoa, &
619 u_c_anorm_drhoanorm_drhoa, u_c_arhoa, u_c_arhoanorm_drhoa, u_c_arhoarhoa, u_c_b, &
620 u_c_b1rhob, u_c_bnorm_drhob, u_c_bnorm_drhobnorm_drhob, u_c_brhob, u_c_brhobnorm_drhob, &
621 u_c_brhobrhob, u_x_a, u_x_a1rhoa, u_x_anorm_drhoa, u_x_anorm_drhoanorm_drhoa, u_x_arhoa, &
622 u_x_arhoanorm_drhoa, u_x_arhoarhoa, u_x_b, u_x_b1rhob, u_x_bnorm_drhob, &
623 u_x_bnorm_drhobnorm_drhob, u_x_brhob, u_x_brhobnorm_drhob, u_x_brhobrhob
624 REAL(kind=dp), DIMENSION(10) :: coeffs
625
626 p_1 = 0.10e1_dp
627 a_1 = 0.31091e-1_dp
628 alpha_1_1 = 0.21370e0_dp
629 beta_1_1 = 0.75957e1_dp
630 beta_2_1 = 0.35876e1_dp
631 beta_3_1 = 0.16382e1_dp
632 beta_4_1 = 0.49294e0_dp
633 p_2 = 0.10e1_dp
634 a_2 = 0.15545e-1_dp
635 alpha_1_2 = 0.20548e0_dp
636 beta_1_2 = 0.141189e2_dp
637 beta_2_2 = 0.61977e1_dp
638 beta_3_2 = 0.33662e1_dp
639 beta_4_2 = 0.62517e0_dp
640 p_3 = 0.10e1_dp
641 a_3 = 0.16887e-1_dp
642 alpha_1_3 = 0.11125e0_dp
643 beta_1_3 = 0.10357e2_dp
644 beta_2_3 = 0.36231e1_dp
645 beta_3_3 = 0.88026e0_dp
646 beta_4_3 = 0.49671e0_dp
647
648 coeffs = b97_coeffs(param)
649 c_x_0 = coeffs(1)
650 c_x_1 = coeffs(2)
651 c_x_2 = coeffs(3)
652 c_cab_0 = coeffs(4)
653 c_cab_1 = coeffs(5)
654 c_cab_2 = coeffs(6)
655 c_css_0 = coeffs(7)
656 c_css_1 = coeffs(8)
657 c_css_2 = coeffs(9)
658
659 scale_c = scale_c_in
660 scale_x = scale_x_in
661 IF (scale_x == -1.0_dp) scale_x = coeffs(10)
662
663 gamma_x = 0.004_dp
664 gamma_c_ss = 0.2_dp
665 gamma_c_ab = 0.006_dp
666
667 SELECT CASE (grad_deriv)
668 CASE default
669 t1 = 3**(0.1e1_dp/0.3e1_dp)
670 t2 = 4**(0.1e1_dp/0.3e1_dp)
671 t3 = t2**2
672 t4 = t1*t3
673 t6 = (0.1e1_dp/pi)**(0.1e1_dp/0.3e1_dp)
674!$OMP DO
675 DO ii = 1, npoints
676 my_rhoa = max(rhoa(ii), 0.0_dp)
677 my_rhob = max(rhob(ii), 0.0_dp)
678 rho = my_rhoa + my_rhob
679 IF (rho > epsilon_rho) THEN
680 my_rhoa = max(my_rhoa, 0.5_dp*epsilon_rho)
681 my_rhob = max(my_rhob, 0.5_dp*epsilon_rho)
682 rho = my_rhoa + my_rhob
683 my_norm_drhoa = norm_drhoa(ii)
684 my_norm_drhob = norm_drhob(ii)
685 t7 = my_rhoa**(0.1e1_dp/0.3e1_dp)
686 t8 = t7*my_rhoa
687 e_lsda_x_a = -0.3e1_dp/0.8e1_dp*t4*t6*t8
688 t12 = 0.1e1_dp/t8
689 s_a = my_norm_drhoa*t12
690 t13 = s_a**2
691 t14 = gamma_x*t13
692 t15 = 0.1e1_dp + t14
693 t16 = 0.1e1_dp/t15
694 u_x_a = t14*t16
695 t18 = c_x_1 + u_x_a*c_x_2
696 gx_a = c_x_0 + u_x_a*t18
697 t20 = my_rhob**(0.1e1_dp/0.3e1_dp)
698 t21 = t20*my_rhob
699 e_lsda_x_b = -0.3e1_dp/0.8e1_dp*t4*t6*t21
700 t25 = 0.1e1_dp/t21
701 s_b = my_norm_drhob*t25
702 t26 = s_b**2
703 t27 = gamma_x*t26
704 t28 = 0.1e1_dp + t27
705 t29 = 0.1e1_dp/t28
706 u_x_b = t27*t29
707 t31 = c_x_1 + u_x_b*c_x_2
708 gx_b = c_x_0 + u_x_b*t31
709 t33 = my_rhoa - my_rhob
710 t34 = 0.1e1_dp/rho
711 chi = t33*t34
712 t35 = 0.1e1_dp/pi
713 t36 = t35*t34
714 t37 = t36**(0.1e1_dp/0.3e1_dp)
715 rs = t4*t37/0.4e1_dp
716 t40 = 0.1e1_dp + alpha_1_1*rs
717 t42 = 0.1e1_dp/a_1
718 t43 = sqrt(rs)
719 t46 = t43*rs
720 t48 = p_1 + 0.1e1_dp
721 t49 = rs**t48
722 t50 = beta_4_1*t49
723 t51 = beta_1_1*t43 + beta_2_1*rs + beta_3_1*t46 + t50
724 t55 = 0.1e1_dp + t42/t51/0.2e1_dp
725 t56 = log(t55)
726 e_c_u_0 = -0.2e1_dp*a_1*t40*t56
727 t60 = 0.1e1_dp + alpha_1_2*rs
728 t62 = 0.1e1_dp/a_2
729 t66 = p_2 + 0.1e1_dp
730 t67 = rs**t66
731 t68 = beta_4_2*t67
732 t69 = beta_1_2*t43 + beta_2_2*rs + beta_3_2*t46 + t68
733 t73 = 0.1e1_dp + t62/t69/0.2e1_dp
734 t74 = log(t73)
735 t78 = 0.1e1_dp + alpha_1_3*rs
736 t80 = 0.1e1_dp/a_3
737 t84 = p_3 + 0.1e1_dp
738 t85 = rs**t84
739 t86 = beta_4_3*t85
740 t87 = beta_1_3*t43 + beta_2_3*rs + beta_3_3*t46 + t86
741 t91 = 0.1e1_dp + t80/t87/0.2e1_dp
742 t92 = log(t91)
743 alpha_c = 0.2e1_dp*a_3*t78*t92
744 t94 = 2**(0.1e1_dp/0.3e1_dp)
745 t97 = 1/(2*t94 - 2)
746 t98 = 0.1e1_dp + chi
747 t99 = t98**(0.1e1_dp/0.3e1_dp)
748 t101 = 0.1e1_dp - chi
749 t102 = t101**(0.1e1_dp/0.3e1_dp)
750 f = (t99*t98 + t102*t101 - 0.2e1_dp)*t97
751 t105 = alpha_c*f
752 t106 = 0.9e1_dp/0.8e1_dp/t97
753 t107 = chi**2
754 t108 = t107**2
755 t110 = t106*(0.1e1_dp - t108)
756 t112 = -0.2e1_dp*a_2*t60*t74 - e_c_u_0
757 t113 = t112*f
758 epsilon_c_unif = e_c_u_0 + t105*t110 + t113*t108
759 t116 = t35/my_rhoa
760 t117 = t116**(0.1e1_dp/0.3e1_dp)
761 rs_a = t4*t117/0.4e1_dp
762 t120 = 0.1e1_dp + alpha_1_2*rs_a
763 t122 = sqrt(rs_a)
764 t125 = t122*rs_a
765 t127 = rs_a**t66
766 t128 = beta_4_2*t127
767 t129 = beta_1_2*t122 + beta_2_2*rs_a + beta_3_2*t125 + t128
768 t133 = 0.1e1_dp + t62/t129/0.2e1_dp
769 t134 = log(t133)
770 epsilon_c_unif_a = -0.2e1_dp*a_2*t120*t134
771 t138 = t35/my_rhob
772 t139 = t138**(0.1e1_dp/0.3e1_dp)
773 rs_b = t4*t139/0.4e1_dp
774 t142 = 0.1e1_dp + alpha_1_2*rs_b
775 t144 = sqrt(rs_b)
776 t147 = t144*rs_b
777 t149 = rs_b**t66
778 t150 = beta_4_2*t149
779 t151 = beta_1_2*t144 + beta_2_2*rs_b + beta_3_2*t147 + t150
780 t155 = 0.1e1_dp + t62/t151/0.2e1_dp
781 t156 = log(t155)
782 epsilon_c_unif_b = -0.2e1_dp*a_2*t142*t156
783 s_a_2 = t13
784 s_b_2 = t26
785 s_avg_2 = s_a_2/0.2e1_dp + s_b_2/0.2e1_dp
786 e_lsda_c_a = epsilon_c_unif_a*my_rhoa
787 e_lsda_c_b = epsilon_c_unif_b*my_rhob
788 t160 = gamma_c_ab*s_avg_2
789 t161 = 0.1e1_dp + t160
790 t162 = 0.1e1_dp/t161
791 u_c_ab = t160*t162
792 t163 = gamma_c_ss*s_a_2
793 t164 = 0.1e1_dp + t163
794 t165 = 0.1e1_dp/t164
795 u_c_a = t163*t165
796 t166 = gamma_c_ss*s_b_2
797 t167 = 0.1e1_dp + t166
798 t168 = 0.1e1_dp/t167
799 u_c_b = t166*t168
800 e_lsda_c_ab = epsilon_c_unif*rho - e_lsda_c_a - e_lsda_c_b
801 t170 = c_cab_1 + u_c_ab*c_cab_2
802 gc_ab = c_cab_0 + u_c_ab*t170
803 t173 = c_css_1 + u_c_a*c_css_2
804 gc_a = c_css_0 + u_c_a*t173
805 t176 = c_css_1 + u_c_b*c_css_2
806 gc_b = c_css_0 + u_c_b*t176
807
808 exc = scale_x*(e_lsda_x_a*gx_a + e_lsda_x_b*gx_b) + scale_c &
809 *(e_lsda_c_ab*gc_ab + e_lsda_c_a*gc_a + e_lsda_c_b*gc_b)
810 e_0(ii) = e_0(ii) + exc
811
812 IF (grad_deriv /= 0) THEN
813
814 e_lsda_x_arhoa = -t4*t6*t7/0.2e1_dp
815 t186 = my_rhoa**2
816 t188 = 0.1e1_dp/t7/t186
817 s_arhoa = -0.4e1_dp/0.3e1_dp*my_norm_drhoa*t188
818 t191 = gamma_x*s_a
819 t192 = t16*s_arhoa
820 t194 = gamma_x**2
821 t196 = t194*s_a_2*s_a
822 t197 = t15**2
823 t198 = 0.1e1_dp/t197
824 t199 = t198*s_arhoa
825 u_x_arhoa = 0.2e1_dp*t191*t192 - 0.2e1_dp*t196*t199
826 gx_arhoa = u_x_arhoa*t18 + u_x_a*u_x_arhoa*c_x_2
827 t207 = rho**2
828 t208 = 0.1e1_dp/t207
829 t209 = t33*t208
830 chirhoa = t34 - t209
831 t210 = t37**2
832 t212 = 0.1e1_dp/t210*t35
833 rsrhoa = -t4*t212*t208/0.12e2_dp
834 t216 = a_1*alpha_1_1
835 t220 = t51**2
836 t221 = 0.1e1_dp/t220
837 t222 = t40*t221
838 t223 = 0.1e1_dp/t43
839 t224 = beta_1_1*t223
840 t228 = beta_3_1*t43
841 t232 = 0.1e1_dp/rs
842 t235 = t224*rsrhoa/0.2e1_dp + beta_2_1*rsrhoa + 0.3e1_dp/ &
843 0.2e1_dp*t228*rsrhoa + t50*t48*rsrhoa*t232
844 t236 = 0.1e1_dp/t55
845 t237 = t235*t236
846 e_c_u_0rhoa = -0.2e1_dp*t216*rsrhoa*t56 + t222*t237
847 t239 = a_2*alpha_1_2
848 t243 = t69**2
849 t244 = 0.1e1_dp/t243
850 t245 = t60*t244
851 t246 = beta_1_2*t223
852 t250 = beta_3_2*t43
853 t256 = t246*rsrhoa/0.2e1_dp + beta_2_2*rsrhoa + 0.3e1_dp/ &
854 0.2e1_dp*t250*rsrhoa + t68*t66*rsrhoa*t232
855 t257 = 0.1e1_dp/t73
856 t258 = t256*t257
857 e_c_u_1rhoa = -0.2e1_dp*t239*rsrhoa*t74 + t245*t258
858 t260 = a_3*alpha_1_3
859 t264 = t87**2
860 t265 = 0.1e1_dp/t264
861 t266 = t78*t265
862 t267 = beta_1_3*t223
863 t271 = beta_3_3*t43
864 t277 = t267*rsrhoa/0.2e1_dp + beta_2_3*rsrhoa + 0.3e1_dp/ &
865 0.2e1_dp*t271*rsrhoa + t86*t84*rsrhoa*t232
866 t278 = 0.1e1_dp/t91
867 t279 = t277*t278
868 alpha_crhoa = 0.2e1_dp*t260*rsrhoa*t92 - t266*t279
869 frhoa = (0.4e1_dp/0.3e1_dp*t99*chirhoa - 0.4e1_dp/0.3e1_dp &
870 *t102*chirhoa)*t97
871 t285 = alpha_crhoa*f
872 t287 = alpha_c*frhoa
873 t289 = t107*chi
874 t290 = t106*t289
875 t291 = t290*chirhoa
876 t293 = 0.4e1_dp*t105*t291
877 t294 = e_c_u_1rhoa - e_c_u_0rhoa
878 t295 = t294*f
879 t297 = t112*frhoa
880 t299 = t289*chirhoa
881 t301 = 0.4e1_dp*t113*t299
882 epsilon_c_unifrhoa = e_c_u_0rhoa + t285*t110 + t287*t110 - &
883 t293 + t295*t108 + t297*t108 + t301
884 t302 = t117**2
885 t304 = 0.1e1_dp/t302*t35
886 rs_arhoa = -t4*t304/t186/0.12e2_dp
887 t312 = t129**2
888 t313 = 0.1e1_dp/t312
889 t314 = t120*t313
890 t315 = 0.1e1_dp/t122
891 t316 = beta_1_2*t315
892 t320 = beta_3_2*t122
893 t324 = 0.1e1_dp/rs_a
894 t327 = t316*rs_arhoa/0.2e1_dp + beta_2_2*rs_arhoa + 0.3e1_dp &
895 /0.2e1_dp*t320*rs_arhoa + t128*t66*rs_arhoa*t324
896 t328 = 0.1e1_dp/t133
897 epsilon_c_unif_arhoa = -0.2e1_dp*t239*rs_arhoa*t134 + t314* &
898 t327*t328
899 s_a_2rhoa = 0.2e1_dp*s_a*s_arhoa
900 s_avg_2rhoa = s_a_2rhoa/0.2e1_dp
901 e_lsda_c_arhoa = epsilon_c_unif_arhoa*my_rhoa + epsilon_c_unif_a
902 t336 = gamma_c_ab**2
903 t337 = t336*s_avg_2
904 t338 = t161**2
905 t339 = 0.1e1_dp/t338
906 u_c_abrhoa = gamma_c_ab*s_avg_2rhoa*t162 - t337*t339*s_avg_2rhoa
907 t344 = gamma_c_ss**2
908 t345 = t344*s_a_2
909 t346 = t164**2
910 t347 = 0.1e1_dp/t346
911 u_c_arhoa = gamma_c_ss*s_a_2rhoa*t165 - t345*t347*s_a_2rhoa
912 e_lsda_c_abrhoa = epsilon_c_unifrhoa*rho + epsilon_c_unif - &
913 e_lsda_c_arhoa
914 gc_abrhoa = u_c_abrhoa*t170 + u_c_ab*u_c_abrhoa*c_cab_2
915 gc_arhoa = u_c_arhoa*t173 + u_c_a*u_c_arhoa*c_css_2
916
917 exc_rhoa = scale_x*(e_lsda_x_arhoa*gx_a + e_lsda_x_a* &
918 gx_arhoa) + scale_c*(e_lsda_c_abrhoa*gc_ab + e_lsda_c_ab* &
919 gc_abrhoa + e_lsda_c_arhoa*gc_a + e_lsda_c_a*gc_arhoa)
920 e_ra(ii) = e_ra(ii) + exc_rhoa
921
922 e_lsda_x_brhob = -t4*t6*t20/0.2e1_dp
923 t365 = my_rhob**2
924 t367 = 0.1e1_dp/t20/t365
925 s_brhob = -0.4e1_dp/0.3e1_dp*my_norm_drhob*t367
926 t370 = gamma_x*s_b
927 t371 = t29*s_brhob
928 t374 = t194*s_b_2*s_b
929 t375 = t28**2
930 t376 = 0.1e1_dp/t375
931 t377 = t376*s_brhob
932 u_x_brhob = 0.2e1_dp*t370*t371 - 0.2e1_dp*t374*t377
933 gx_brhob = u_x_brhob*t31 + u_x_b*u_x_brhob*c_x_2
934 chirhob = -t34 - t209
935 rsrhob = rsrhoa
936 t396 = t224*rsrhob/0.2e1_dp + beta_2_1*rsrhob + 0.3e1_dp/ &
937 0.2e1_dp*t228*rsrhob + t50*t48*rsrhob*t232
938 e_c_u_0rhob = -0.2e1_dp*t216*rsrhob*t56 + t222*t396*t236
939 t410 = t246*rsrhob/0.2e1_dp + beta_2_2*rsrhob + 0.3e1_dp/ &
940 0.2e1_dp*t250*rsrhob + t68*t66*rsrhob*t232
941 e_c_u_1rhob = -0.2e1_dp*t239*rsrhob*t74 + t245*t410*t257
942 t424 = t267*rsrhob/0.2e1_dp + beta_2_3*rsrhob + 0.3e1_dp/ &
943 0.2e1_dp*t271*rsrhob + t86*t84*rsrhob*t232
944 alpha_crhob = 0.2e1_dp*t260*rsrhob*t92 - t266*t424*t278
945 frhob = (0.4e1_dp/0.3e1_dp*t99*chirhob - 0.4e1_dp/0.3e1_dp &
946 *t102*chirhob)*t97
947 t431 = alpha_crhob*f
948 t433 = alpha_c*frhob
949 t435 = t290*chirhob
950 t437 = 0.4e1_dp*t105*t435
951 t438 = e_c_u_1rhob - e_c_u_0rhob
952 t439 = t438*f
953 t441 = t112*frhob
954 t443 = t289*chirhob
955 t445 = 0.4e1_dp*t113*t443
956 epsilon_c_unifrhob = e_c_u_0rhob + t431*t110 + t433*t110 - &
957 t437 + t439*t108 + t441*t108 + t445
958 t446 = t139**2
959 t448 = 0.1e1_dp/t446*t35
960 rs_brhob = -t4*t448/t365/0.12e2_dp
961 t456 = t151**2
962 t457 = 0.1e1_dp/t456
963 t458 = t142*t457
964 t459 = 0.1e1_dp/t144
965 t460 = beta_1_2*t459
966 t464 = beta_3_2*t144
967 t468 = 0.1e1_dp/rs_b
968 t471 = t460*rs_brhob/0.2e1_dp + beta_2_2*rs_brhob + 0.3e1_dp &
969 /0.2e1_dp*rs_brhob*t464 + t150*t66*rs_brhob*t468
970 t472 = 0.1e1_dp/t155
971 epsilon_c_unif_brhob = -0.2e1_dp*t239*rs_brhob*t156 + t458* &
972 t471*t472
973 s_b_2rhob = 0.2e1_dp*s_b*s_brhob
974 s_avg_2rhob = s_b_2rhob/0.2e1_dp
975 e_lsda_c_brhob = epsilon_c_unif_brhob*my_rhob + epsilon_c_unif_b
976 t480 = t339*s_avg_2rhob
977 u_c_abrhob = gamma_c_ab*s_avg_2rhob*t162 - t337*t480
978 t484 = t344*s_b_2
979 t485 = t167**2
980 t486 = 0.1e1_dp/t485
981 u_c_brhob = gamma_c_ss*s_b_2rhob*t168 - t484*t486*s_b_2rhob
982 e_lsda_c_abrhob = epsilon_c_unifrhob*rho + epsilon_c_unif - &
983 e_lsda_c_brhob
984 gc_abrhob = u_c_abrhob*t170 + u_c_ab*u_c_abrhob*c_cab_2
985 gc_brhob = u_c_brhob*t176 + u_c_b*u_c_brhob*c_css_2
986
987 exc_rhob = scale_x*(e_lsda_x_brhob*gx_b + e_lsda_x_b* &
988 gx_brhob) + scale_c*(e_lsda_c_abrhob*gc_ab + e_lsda_c_ab* &
989 gc_abrhob + e_lsda_c_brhob*gc_b + e_lsda_c_b*gc_brhob)
990 e_rb(ii) = e_rb(ii) + exc_rhob
991
992 s_anorm_drhoa = t12
993 u_x_anorm_drhoa = 0.2e1_dp*t191*t16*s_anorm_drhoa - 0.2e1_dp &
994 *t196*t198*s_anorm_drhoa
995 gx_anorm_drhoa = u_x_anorm_drhoa*t18 + u_x_a*u_x_anorm_drhoa*c_x_2
996 s_a_2norm_drhoa = 0.2e1_dp*s_a*s_anorm_drhoa
997 s_avg_2norm_drhoa = s_a_2norm_drhoa/0.2e1_dp
998 t512 = t339*s_avg_2norm_drhoa
999 u_c_abnorm_drhoa = gamma_c_ab*s_avg_2norm_drhoa*t162 - t337*t512
1000 t516 = t347*s_a_2norm_drhoa
1001 u_c_anorm_drhoa = gamma_c_ss*s_a_2norm_drhoa*t165 - t345*t516
1002 gc_abnorm_drhoa = u_c_abnorm_drhoa*t170 + u_c_ab* &
1003 u_c_abnorm_drhoa*c_cab_2
1004 gc_anorm_drhoa = u_c_anorm_drhoa*t173 + u_c_a*u_c_anorm_drhoa &
1005 *c_css_2
1006
1007 exc_norm_drhoa = scale_x*e_lsda_x_a*gx_anorm_drhoa + scale_c* &
1008 (e_lsda_c_ab*gc_abnorm_drhoa + e_lsda_c_a*gc_anorm_drhoa)
1009 e_ndra(ii) = e_ndra(ii) + exc_norm_drhoa
1010
1011 s_bnorm_drhob = t25
1012 u_x_bnorm_drhob = 0.2e1_dp*t370*t29*s_bnorm_drhob - 0.2e1_dp &
1013 *t374*t376*s_bnorm_drhob
1014 gx_bnorm_drhob = u_x_bnorm_drhob*t31 + u_x_b*u_x_bnorm_drhob*c_x_2
1015 s_b_2norm_drhob = 0.2e1_dp*s_b*s_bnorm_drhob
1016 s_avg_2norm_drhob = s_b_2norm_drhob/0.2e1_dp
1017 t539 = t339*s_avg_2norm_drhob
1018 u_c_abnorm_drhob = gamma_c_ab*s_avg_2norm_drhob*t162 - t337*t539
1019 t543 = t486*s_b_2norm_drhob
1020 u_c_bnorm_drhob = gamma_c_ss*s_b_2norm_drhob*t168 - t484*t543
1021 gc_abnorm_drhob = u_c_abnorm_drhob*t170 + u_c_ab* &
1022 u_c_abnorm_drhob*c_cab_2
1023 gc_bnorm_drhob = u_c_bnorm_drhob*t176 + u_c_b*u_c_bnorm_drhob &
1024 *c_css_2
1025
1026 exc_norm_drhob = scale_x*e_lsda_x_b*gx_bnorm_drhob + scale_c* &
1027 (e_lsda_c_ab*gc_abnorm_drhob + e_lsda_c_b*gc_bnorm_drhob)
1028 e_ndrb(ii) = e_ndrb(ii) + exc_norm_drhob
1029
1030 IF (grad_deriv > 1) THEN
1031 t555 = t7**2
1032 t560 = t186*my_rhoa
1033 s_arhoarhoa = 0.28e2_dp/0.9e1_dp*my_norm_drhoa/t7/t560
1034 t564 = s_arhoa**2
1035 t568 = t194*s_a_2
1036 t575 = t194*gamma_x
1037 t576 = s_a_2**2
1038 t577 = t575*t576
1039 t579 = 0.1e1_dp/t197/t15
1040 u_x_arhoarhoa = 0.2e1_dp*gamma_x*t564*t16 - 0.10e2_dp*t568 &
1041 *t198*t564 + 0.2e1_dp*t191*t16*s_arhoarhoa + 0.8e1_dp* &
1042 t577*t579*t564 - 0.2e1_dp*t196*t198*s_arhoarhoa
1043 u_x_a1rhoa = u_x_arhoa
1044 t600 = 0.1e1_dp/t207/rho
1045 t601 = t33*t600
1046 chirhoarhoa = -0.2e1_dp*t208 + 0.2e1_dp*t601
1047 t605 = 0.3141592654e1_dp**2
1048 t606 = 0.1e1_dp/t605
1049 t608 = t207**2
1050 rsrhoarhoa = -t4/t210/t36*t606/t608/0.18e2_dp + &
1051 t4*t212*t600/0.6e1_dp
1052 t619 = alpha_1_1*rsrhoa
1053 t621 = t221*t235*t236
1054 t626 = t40/t220/t51
1055 t627 = t235**2
1056 t631 = 0.1e1_dp/t46
1057 t632 = beta_1_1*t631
1058 t633 = rsrhoa**2
1059 t639 = beta_3_1*t223
1060 t644 = t48**2
1061 t646 = rs**2
1062 t647 = 0.1e1_dp/t646
1063 t659 = t220**2
1064 t661 = t40/t659
1065 t662 = t55**2
1066 t663 = 0.1e1_dp/t662
1067 e_c_u_0rhoarhoa = -0.2e1_dp*t216*rsrhoarhoa*t56 + 0.2e1_dp* &
1068 t619*t621 - 0.2e1_dp*t626*t627*t236 + t222*(-t632*t633 &
1069 /0.4e1_dp + t224*rsrhoarhoa/0.2e1_dp + beta_2_1*rsrhoarhoa + &
1070 0.3e1_dp/0.4e1_dp*t639*t633 + 0.3e1_dp/0.2e1_dp*t228* &
1071 rsrhoarhoa + t50*t644*t633*t647 + t50*t48*rsrhoarhoa* &
1072 t232 - t50*t48*t633*t647)*t236 + t661*t627*t663*t42/ &
1073 0.2e1_dp
1074 e_c_u_01rhoa = e_c_u_0rhoa
1075 t671 = alpha_1_2*rsrhoa
1076 t673 = t244*t256*t257
1077 t678 = t60/t243/t69
1078 t679 = t256**2
1079 t683 = beta_1_2*t631
1080 t689 = beta_3_2*t223
1081 t694 = t66**2
1082 t707 = t243**2
1083 t709 = t60/t707
1084 t710 = t73**2
1085 t711 = 0.1e1_dp/t710
1086 t719 = alpha_1_3*rsrhoa
1087 t721 = t265*t277*t278
1088 t726 = t78/t264/t87
1089 t727 = t277**2
1090 t731 = beta_1_3*t631
1091 t737 = beta_3_3*t223
1092 t742 = t84**2
1093 t755 = t264**2
1094 t757 = t78/t755
1095 t758 = t91**2
1096 t759 = 0.1e1_dp/t758
1097 alpha_c1rhoa = alpha_crhoa
1098 t764 = t99**2
1099 t765 = 0.1e1_dp/t764
1100 t766 = chirhoa**2
1101 t771 = t102**2
1102 t772 = 0.1e1_dp/t771
1103 frhoarhoa = (0.4e1_dp/0.9e1_dp*t765*t766 + 0.4e1_dp/ &
1104 0.3e1_dp*t99*chirhoarhoa + 0.4e1_dp/0.9e1_dp*t772*t766 - &
1105 0.4e1_dp/0.3e1_dp*t102*chirhoarhoa)*t97
1106 f1rhoa = frhoa
1107 t790 = alpha_c1rhoa*f
1108 t793 = alpha_c*f1rhoa
1109 t796 = t106*t107
1110 t811 = e_c_u_1rhoa - e_c_u_01rhoa
1111 t818 = t811*f
1112 t821 = t112*f1rhoa
1113 t830 = -0.4e1_dp*t105*t290*chirhoarhoa + (-0.2e1_dp*t239* &
1114 rsrhoarhoa*t74 + 0.2e1_dp*t671*t673 - 0.2e1_dp*t678*t679 &
1115 *t257 + t245*(-t683*t633/0.4e1_dp + t246*rsrhoarhoa/ &
1116 0.2e1_dp + beta_2_2*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t689*t633 &
1117 + 0.3e1_dp/0.2e1_dp*t250*rsrhoarhoa + t68*t694*t633* &
1118 t647 + t68*t66*rsrhoarhoa*t232 - t68*t66*t633*t647)* &
1119 t257 + t709*t679*t711*t62/0.2e1_dp - e_c_u_0rhoarhoa)*f* &
1120 t108 + t294*f1rhoa*t108 + 0.4e1_dp*t295*t299 + t811*frhoa &
1121 *t108 + t112*frhoarhoa*t108 + 0.4e1_dp*t297*t299 + 0.4e1_dp &
1122 *t818*t299 + 0.4e1_dp*t821*t299 + 0.12e2_dp*t113*t107* &
1123 t766 + 0.4e1_dp*t113*t289*chirhoarhoa
1124 epsilon_c_unif1rhoa = e_c_u_01rhoa + t790*t110 + t793*t110 - &
1125 t293 + t818*t108 + t821*t108 + t301
1126 t838 = t186**2
1127 rs_arhoarhoa = -t4/t302/t116*t606/t838/0.18e2_dp + &
1128 t4*t304/t560/0.6e1_dp
1129 t858 = t327**2
1130 t864 = rs_arhoa**2
1131 t876 = rs_a**2
1132 t877 = 0.1e1_dp/t876
1133 t889 = t312**2
1134 t892 = t133**2
1135 epsilon_c_unif_a1rhoa = epsilon_c_unif_arhoa
1136 s_a_2rhoarhoa = 0.2e1_dp*t564 + 0.2e1_dp*s_a*s_arhoarhoa
1137 s_a_21rhoa = s_a_2rhoa
1138 s_avg_2rhoarhoa = s_a_2rhoarhoa/0.2e1_dp
1139 s_avg_21rhoa = s_a_21rhoa/0.2e1_dp
1140 e_lsda_c_arhoarhoa = (-0.2e1_dp*t239*rs_arhoarhoa*t134 + &
1141 0.2e1_dp*alpha_1_2*rs_arhoa*t313*t327*t328 - 0.2e1_dp* &
1142 t120/t312/t129*t858*t328 + t314*(-beta_1_2/t125*t864/ &
1143 0.4e1_dp + t316*rs_arhoarhoa/0.2e1_dp + beta_2_2*rs_arhoarhoa &
1144 + 0.3e1_dp/0.4e1_dp*beta_3_2*t315*t864 + 0.3e1_dp/ &
1145 0.2e1_dp*t320*rs_arhoarhoa + t128*t694*t864*t877 + t128* &
1146 t66*rs_arhoarhoa*t324 - t128*t66*t864*t877)*t328 + t120 &
1147 /t889*t858/t892*t62/0.2e1_dp)*my_rhoa + epsilon_c_unif_arhoa &
1148 + epsilon_c_unif_a1rhoa
1149 e_lsda_c_a1rhoa = epsilon_c_unif_a1rhoa*my_rhoa + epsilon_c_unif_a
1150 t906 = t336*s_avg_2rhoa
1151 t907 = t339*s_avg_21rhoa
1152 t911 = t336*gamma_c_ab*s_avg_2
1153 t913 = 0.1e1_dp/t338/t161
1154 t914 = t913*s_avg_2rhoa
1155 u_c_abrhoarhoa = gamma_c_ab*s_avg_2rhoarhoa*t162 - 0.2e1_dp* &
1156 t906*t907 + 0.2e1_dp*t911*t914*s_avg_21rhoa - t337*t339* &
1157 s_avg_2rhoarhoa
1158 u_c_ab1rhoa = gamma_c_ab*s_avg_21rhoa*t162 - t337*t907
1159 t925 = t344*s_a_2rhoa
1160 t926 = t347*s_a_21rhoa
1161 t929 = t344*gamma_c_ss
1162 t930 = t929*s_a_2
1163 t932 = 0.1e1_dp/t346/t164
1164 t933 = t932*s_a_2rhoa
1165 u_c_arhoarhoa = gamma_c_ss*s_a_2rhoarhoa*t165 - 0.2e1_dp* &
1166 t925*t926 + 0.2e1_dp*t930*t933*s_a_21rhoa - t345*t347* &
1167 s_a_2rhoarhoa
1168 u_c_a1rhoa = gamma_c_ss*s_a_21rhoa*t165 - t345*t926
1169
1170 exc_rhoa_rhoa = scale_x*(-t4*t6/t555*gx_a/0.6e1_dp &
1171 + e_lsda_x_arhoa*(u_x_a1rhoa*t18 + u_x_a*u_x_a1rhoa*c_x_2) &
1172 + e_lsda_x_arhoa*gx_arhoa + e_lsda_x_a*(u_x_arhoarhoa*t18 + &
1173 0.2e1_dp*u_x_arhoa*u_x_a1rhoa*c_x_2 + u_x_a*u_x_arhoarhoa* &
1174 c_x_2)) + scale_c*(((e_c_u_0rhoarhoa + (0.2e1_dp*t260* &
1175 rsrhoarhoa*t92 - 0.2e1_dp*t719*t721 + 0.2e1_dp*t726*t727* &
1176 t278 - t266*(-t731*t633/0.4e1_dp + t267*rsrhoarhoa/ &
1177 0.2e1_dp + beta_2_3*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t737*t633 &
1178 + 0.3e1_dp/0.2e1_dp*t271*rsrhoarhoa + t86*t742*t633* &
1179 t647 + t86*t84*rsrhoarhoa*t232 - t86*t84*t633*t647)* &
1180 t278 - t757*t727*t759*t80/0.2e1_dp)*f*t110 + alpha_crhoa &
1181 *f1rhoa*t110 - 0.4e1_dp*t285*t291 + alpha_c1rhoa*frhoa* &
1182 t110 + alpha_c*frhoarhoa*t110 - 0.4e1_dp*t287*t291 - &
1183 0.4e1_dp*t790*t291 - 0.4e1_dp*t793*t291 - 0.12e2_dp*t105* &
1184 t796*t766 + t830)*rho + epsilon_c_unifrhoa + &
1185 epsilon_c_unif1rhoa - e_lsda_c_arhoarhoa)*gc_ab + e_lsda_c_abrhoa &
1186 *(u_c_ab1rhoa*t170 + u_c_ab*u_c_ab1rhoa*c_cab_2) + ( &
1187 epsilon_c_unif1rhoa*rho + epsilon_c_unif - e_lsda_c_a1rhoa)* &
1188 gc_abrhoa + e_lsda_c_ab*(u_c_abrhoarhoa*t170 + 0.2e1_dp* &
1189 u_c_abrhoa*u_c_ab1rhoa*c_cab_2 + u_c_ab*u_c_abrhoarhoa* &
1190 c_cab_2) + e_lsda_c_arhoarhoa*gc_a + e_lsda_c_arhoa*(u_c_a1rhoa &
1191 *t173 + u_c_a*u_c_a1rhoa*c_css_2) + e_lsda_c_a1rhoa*gc_arhoa &
1192 + e_lsda_c_a*(u_c_arhoarhoa*t173 + 0.2e1_dp*u_c_arhoa* &
1193 u_c_a1rhoa*c_css_2 + u_c_a*u_c_arhoarhoa*c_css_2))
1194 e_ra_ra(ii) = e_ra_ra(ii) + exc_rhoa_rhoa
1195
1196 chirhoarhob = 0.2e1_dp*t601
1197 rsrhoarhob = rsrhoarhoa
1198 t974 = t221*t396*t236
1199 t976 = alpha_1_1*rsrhob
1200 t981 = rsrhoa*rsrhob
1201 t993 = rsrhob*t647*rsrhoa
1202 e_c_u_0rhoarhob = -0.2e1_dp*t216*rsrhoarhob*t56 + t619* &
1203 t974 + t976*t621 - 0.2e1_dp*t626*t237*t396 + t222*(-t632* &
1204 t981/0.4e1_dp + t224*rsrhoarhob/0.2e1_dp + beta_2_1* &
1205 rsrhoarhob + 0.3e1_dp/0.4e1_dp*t639*t981 + 0.3e1_dp/0.2e1_dp &
1206 *t228*rsrhoarhob + t50*t644*t993 + t50*t48*rsrhoarhob* &
1207 t232 - t50*t48*t993)*t236 + t661*t235*t663*t42*t396/ &
1208 0.2e1_dp
1209 t1012 = t244*t410*t257
1210 t1014 = alpha_1_2*rsrhob
1211 t1047 = t265*t424*t278
1212 t1049 = alpha_1_3*rsrhob
1213 frhoarhob = (0.4e1_dp/0.9e1_dp*t765*chirhoa*chirhob + &
1214 0.4e1_dp/0.3e1_dp*t99*chirhoarhob + 0.4e1_dp/0.9e1_dp*t772 &
1215 *chirhoa*chirhob - 0.4e1_dp/0.3e1_dp*t102*chirhoarhob)* &
1216 t97
1217 t1107 = t107*chirhoa*chirhob
1218 t1136 = -0.4e1_dp*t105*t290*chirhoarhob + (-0.2e1_dp*t239 &
1219 *rsrhoarhob*t74 + t671*t1012 + t1014*t673 - 0.2e1_dp*t678* &
1220 t258*t410 + t245*(-t683*t981/0.4e1_dp + t246*rsrhoarhob/ &
1221 0.2e1_dp + beta_2_2*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t689* &
1222 t981 + 0.3e1_dp/0.2e1_dp*t250*rsrhoarhob + t68*t694*t993 + &
1223 t68*t66*rsrhoarhob*t232 - t68*t66*t993)*t257 + t709* &
1224 t256*t711*t62*t410/0.2e1_dp - e_c_u_0rhoarhob)*f*t108 + &
1225 t294*frhob*t108 + 0.4e1_dp*t295*t443 + t438*frhoa*t108 + &
1226 t112*frhoarhob*t108 + 0.4e1_dp*t297*t443 + 0.4e1_dp*t439 &
1227 *t299 + 0.4e1_dp*t441*t299 + 0.12e2_dp*t113*t1107 + &
1228 0.4e1_dp*t113*t289*chirhoarhob
1229 u_c_abrhoarhob = -0.2e1_dp*t906*t480 + 0.2e1_dp*t911*t914 &
1230 *s_avg_2rhob
1231
1232 exc_rhoa_rhob = scale_c*(((e_c_u_0rhoarhob + (0.2e1_dp*t260* &
1233 rsrhoarhob*t92 - t719*t1047 - t1049*t721 + 0.2e1_dp*t726* &
1234 t279*t424 - t266*(-t731*t981/0.4e1_dp + t267*rsrhoarhob/ &
1235 0.2e1_dp + beta_2_3*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t737*t981 &
1236 + 0.3e1_dp/0.2e1_dp*t271*rsrhoarhob + t86*t742*t993 + t86 &
1237 *t84*rsrhoarhob*t232 - t86*t84*t993)*t278 - t757*t277 &
1238 *t759*t80*t424/0.2e1_dp)*f*t110 + alpha_crhoa*frhob* &
1239 t110 - 0.4e1_dp*t285*t435 + alpha_crhob*frhoa*t110 + alpha_c &
1240 *frhoarhob*t110 - 0.4e1_dp*t287*t435 - 0.4e1_dp*t431* &
1241 t291 - 0.4e1_dp*t433*t291 - 0.12e2_dp*t105*t106*t1107 + &
1242 t1136)*rho + epsilon_c_unifrhoa + epsilon_c_unifrhob)*gc_ab + &
1243 e_lsda_c_abrhoa*gc_abrhob + e_lsda_c_abrhob*gc_abrhoa + &
1244 e_lsda_c_ab*(u_c_abrhoarhob*t170 + 0.2e1_dp*u_c_abrhoa* &
1245 u_c_abrhob*c_cab_2 + u_c_ab*u_c_abrhoarhob*c_cab_2))
1246 e_ra_rb(ii) = e_ra_rb(ii) + exc_rhoa_rhob
1247
1248 t1152 = t20**2
1249 t1157 = t365*my_rhob
1250 s_brhobrhob = 0.28e2_dp/0.9e1_dp*my_norm_drhob/t20/t1157
1251 t1161 = s_brhob**2
1252 t1165 = t194*s_b_2
1253 t1172 = s_b_2**2
1254 t1173 = t575*t1172
1255 t1175 = 0.1e1_dp/t375/t28
1256 u_x_brhobrhob = 0.2e1_dp*gamma_x*t1161*t29 - 0.10e2_dp* &
1257 t1165*t376*t1161 + 0.2e1_dp*t370*t29*s_brhobrhob + &
1258 0.8e1_dp*t1173*t1175*t1161 - 0.2e1_dp*t374*t376* &
1259 s_brhobrhob
1260 u_x_b1rhob = u_x_brhob
1261 chirhobrhob = 0.2e1_dp*t208 + 0.2e1_dp*t601
1262 rsrhobrhob = rsrhoarhob
1263 t1201 = t396**2
1264 t1205 = rsrhob**2
1265 e_c_u_0rhobrhob = -0.2e1_dp*t216*rsrhobrhob*t56 + 0.2e1_dp* &
1266 t976*t974 - 0.2e1_dp*t626*t1201*t236 + t222*(-t632* &
1267 t1205/0.4e1_dp + t224*rsrhobrhob/0.2e1_dp + beta_2_1* &
1268 rsrhobrhob + 0.3e1_dp/0.4e1_dp*t639*t1205 + 0.3e1_dp/ &
1269 0.2e1_dp*t228*rsrhobrhob + t50*t644*t1205*t647 + t50*t48 &
1270 *rsrhobrhob*t232 - t50*t48*t1205*t647)*t236 + t661* &
1271 t1201*t663*t42/0.2e1_dp
1272 e_c_u_01rhob = e_c_u_0rhob
1273 t1236 = t410**2
1274 t1270 = t424**2
1275 alpha_c1rhob = alpha_crhob
1276 t1299 = chirhob**2
1277 frhobrhob = (0.4e1_dp/0.9e1_dp*t765*t1299 + 0.4e1_dp/ &
1278 0.3e1_dp*t99*chirhobrhob + 0.4e1_dp/0.9e1_dp*t772*t1299 - &
1279 0.4e1_dp/0.3e1_dp*t102*chirhobrhob)*t97
1280 f1rhob = frhob
1281 t1321 = alpha_c1rhob*f
1282 t1324 = alpha_c*f1rhob
1283 t1341 = e_c_u_1rhob - e_c_u_01rhob
1284 t1348 = t1341*f
1285 t1351 = t112*f1rhob
1286 t1360 = -0.4e1_dp*t105*t290*chirhobrhob + (-0.2e1_dp*t239 &
1287 *rsrhobrhob*t74 + 0.2e1_dp*t1014*t1012 - 0.2e1_dp*t678* &
1288 t1236*t257 + t245*(-t683*t1205/0.4e1_dp + t246*rsrhobrhob &
1289 /0.2e1_dp + beta_2_2*rsrhobrhob + 0.3e1_dp/0.4e1_dp*t689* &
1290 t1205 + 0.3e1_dp/0.2e1_dp*t250*rsrhobrhob + t68*t694*t1205 &
1291 *t647 + t68*t66*rsrhobrhob*t232 - t68*t66*t1205*t647) &
1292 *t257 + t709*t1236*t711*t62/0.2e1_dp - e_c_u_0rhobrhob)*f &
1293 *t108 + t438*f1rhob*t108 + 0.4e1_dp*t439*t443 + t1341* &
1294 frhob*t108 + t112*frhobrhob*t108 + 0.4e1_dp*t441*t443 + &
1295 0.4e1_dp*t1348*t443 + 0.4e1_dp*t1351*t443 + 0.12e2_dp*t113 &
1296 *t107*t1299 + 0.4e1_dp*t113*t289*chirhobrhob
1297 epsilon_c_unif1rhob = e_c_u_01rhob + t1321*t110 + t1324*t110 - &
1298 t437 + t1348*t108 + t1351*t108 + t445
1299 t1368 = t365**2
1300 rs_brhobrhob = -t4/t446/t138*t606/t1368/0.18e2_dp &
1301 + t4*t448/t1157/0.6e1_dp
1302 t1388 = t471**2
1303 t1394 = rs_brhob**2
1304 t1406 = rs_b**2
1305 t1407 = 0.1e1_dp/t1406
1306 t1419 = t456**2
1307 t1422 = t155**2
1308 epsilon_c_unif_b1rhob = epsilon_c_unif_brhob
1309 s_b_2rhobrhob = 0.2e1_dp*t1161 + 0.2e1_dp*s_b*s_brhobrhob
1310 s_b_21rhob = s_b_2rhob
1311 s_avg_2rhobrhob = s_b_2rhobrhob/0.2e1_dp
1312 s_avg_21rhob = s_b_21rhob/0.2e1_dp
1313 e_lsda_c_brhobrhob = (-0.2e1_dp*t239*rs_brhobrhob*t156 + &
1314 0.2e1_dp*alpha_1_2*rs_brhob*t457*t471*t472 - 0.2e1_dp* &
1315 t142/t456/t151*t1388*t472 + t458*(-beta_1_2/t147*t1394 &
1316 /0.4e1_dp + t460*rs_brhobrhob/0.2e1_dp + beta_2_2* &
1317 rs_brhobrhob + 0.3e1_dp/0.4e1_dp*beta_3_2*t459*t1394 + &
1318 0.3e1_dp/0.2e1_dp*t464*rs_brhobrhob + t150*t694*t1394* &
1319 t1407 + t150*t66*rs_brhobrhob*t468 - t150*t66*t1394* &
1320 t1407)*t472 + t142/t1419*t1388/t1422*t62/0.2e1_dp)* &
1321 my_rhob + epsilon_c_unif_brhob + epsilon_c_unif_b1rhob
1322 e_lsda_c_b1rhob = epsilon_c_unif_b1rhob*my_rhob + epsilon_c_unif_b
1323 t1436 = t336*s_avg_2rhob
1324 t1437 = t339*s_avg_21rhob
1325 t1440 = t913*s_avg_2rhob
1326 u_c_abrhobrhob = gamma_c_ab*s_avg_2rhobrhob*t162 - 0.2e1_dp* &
1327 t1436*t1437 + 0.2e1_dp*t911*t1440*s_avg_21rhob - t337*t339 &
1328 *s_avg_2rhobrhob
1329 u_c_ab1rhob = gamma_c_ab*s_avg_21rhob*t162 - t337*t1437
1330 t1451 = t344*s_b_2rhob
1331 t1452 = t486*s_b_21rhob
1332 t1455 = t929*s_b_2
1333 t1457 = 0.1e1_dp/t485/t167
1334 t1458 = t1457*s_b_2rhob
1335 u_c_brhobrhob = gamma_c_ss*s_b_2rhobrhob*t168 - 0.2e1_dp* &
1336 t1451*t1452 + 0.2e1_dp*t1455*t1458*s_b_21rhob - t484*t486 &
1337 *s_b_2rhobrhob
1338 u_c_b1rhob = gamma_c_ss*s_b_21rhob*t168 - t484*t1452
1339
1340 exc_rhob_rhob = scale_x*(-t4*t6/t1152*gx_b/ &
1341 0.6e1_dp + e_lsda_x_brhob*(u_x_b1rhob*t31 + u_x_b*u_x_b1rhob* &
1342 c_x_2) + e_lsda_x_brhob*gx_brhob + e_lsda_x_b*(u_x_brhobrhob* &
1343 t31 + 0.2e1_dp*u_x_brhob*u_x_b1rhob*c_x_2 + u_x_b* &
1344 u_x_brhobrhob*c_x_2)) + scale_c*(((e_c_u_0rhobrhob + (0.2e1_dp* &
1345 t260*rsrhobrhob*t92 - 0.2e1_dp*t1049*t1047 + 0.2e1_dp* &
1346 t726*t1270*t278 - t266*(-t731*t1205/0.4e1_dp + t267* &
1347 rsrhobrhob/0.2e1_dp + beta_2_3*rsrhobrhob + 0.3e1_dp/0.4e1_dp &
1348 *t737*t1205 + 0.3e1_dp/0.2e1_dp*t271*rsrhobrhob + t86* &
1349 t742*t1205*t647 + t86*t84*rsrhobrhob*t232 - t86*t84* &
1350 t1205*t647)*t278 - t757*t1270*t759*t80/0.2e1_dp)*f* &
1351 t110 + alpha_crhob*f1rhob*t110 - 0.4e1_dp*t431*t435 + &
1352 alpha_c1rhob*frhob*t110 + alpha_c*frhobrhob*t110 - 0.4e1_dp &
1353 *t433*t435 - 0.4e1_dp*t1321*t435 - 0.4e1_dp*t1324*t435 - &
1354 0.12e2_dp*t105*t796*t1299 + t1360)*rho + epsilon_c_unifrhob &
1355 + epsilon_c_unif1rhob - e_lsda_c_brhobrhob)*gc_ab + &
1356 e_lsda_c_abrhob*(u_c_ab1rhob*t170 + u_c_ab*u_c_ab1rhob* &
1357 c_cab_2) + (epsilon_c_unif1rhob*rho + epsilon_c_unif - &
1358 e_lsda_c_b1rhob)*gc_abrhob + e_lsda_c_ab*(u_c_abrhobrhob*t170 &
1359 + 0.2e1_dp*u_c_abrhob*u_c_ab1rhob*c_cab_2 + u_c_ab* &
1360 u_c_abrhobrhob*c_cab_2) + e_lsda_c_brhobrhob*gc_b + &
1361 e_lsda_c_brhob*(u_c_b1rhob*t176 + u_c_b*u_c_b1rhob*c_css_2) &
1362 + e_lsda_c_b1rhob*gc_brhob + e_lsda_c_b*(u_c_brhobrhob*t176 + &
1363 0.2e1_dp*u_c_brhob*u_c_b1rhob*c_css_2 + u_c_b*u_c_brhobrhob &
1364 *c_css_2))
1365 e_rb_rb(ii) = e_rb_rb(ii) + exc_rhob_rhob
1366
1367 s_arhoanorm_drhoa = -0.4e1_dp/0.3e1_dp*t188
1368 u_x_arhoanorm_drhoa = 0.2e1_dp*gamma_x*s_anorm_drhoa*t192 - &
1369 0.10e2_dp*t568*t199*s_anorm_drhoa + 0.2e1_dp*t191*t16* &
1370 s_arhoanorm_drhoa + 0.8e1_dp*t577*t579*s_arhoa*s_anorm_drhoa &
1371 - 0.2e1_dp*t196*t198*s_arhoanorm_drhoa
1372 s_a_2rhoanorm_drhoa = 0.2e1_dp*s_anorm_drhoa*s_arhoa + &
1373 0.2e1_dp*s_a*s_arhoanorm_drhoa
1374 s_avg_2rhoanorm_drhoa = s_a_2rhoanorm_drhoa/0.2e1_dp
1375 u_c_abrhoanorm_drhoa = gamma_c_ab*s_avg_2rhoanorm_drhoa*t162 - &
1376 0.2e1_dp*t906*t512 + 0.2e1_dp*t911*t914*s_avg_2norm_drhoa &
1377 - t337*t339*s_avg_2rhoanorm_drhoa
1378 u_c_arhoanorm_drhoa = gamma_c_ss*s_a_2rhoanorm_drhoa*t165 - &
1379 0.2e1_dp*t925*t516 + 0.2e1_dp*t930*t933*s_a_2norm_drhoa - &
1380 t345*t347*s_a_2rhoanorm_drhoa
1381
1382 exc_rhoa_norm_drhoa = scale_x*(e_lsda_x_arhoa*gx_anorm_drhoa + &
1383 e_lsda_x_a*(u_x_arhoanorm_drhoa*t18 + 0.2e1_dp*u_x_arhoa* &
1384 u_x_anorm_drhoa*c_x_2 + u_x_a*u_x_arhoanorm_drhoa*c_x_2)) + &
1385 scale_c*(e_lsda_c_abrhoa*gc_abnorm_drhoa + e_lsda_c_ab*( &
1386 u_c_abrhoanorm_drhoa*t170 + 0.2e1_dp*u_c_abrhoa* &
1387 u_c_abnorm_drhoa*c_cab_2 + u_c_ab*u_c_abrhoanorm_drhoa*c_cab_2 &
1388 ) + e_lsda_c_arhoa*gc_anorm_drhoa + e_lsda_c_a*( &
1389 u_c_arhoanorm_drhoa*t173 + 0.2e1_dp*u_c_arhoa*u_c_anorm_drhoa &
1390 *c_css_2 + u_c_a*u_c_arhoanorm_drhoa*c_css_2))
1391 e_ra_ndra(ii) = e_ra_ndra(ii) + exc_rhoa_norm_drhoa
1392
1393 u_c_abrhobnorm_drhoa = -0.2e1_dp*t1436*t512 + 0.2e1_dp*t911 &
1394 *t1440*s_avg_2norm_drhoa
1395
1396 exc_rhob_norm_drhoa = scale_c*(e_lsda_c_abrhob*gc_abnorm_drhoa &
1397 + e_lsda_c_ab*(u_c_abrhobnorm_drhoa*t170 + 0.2e1_dp* &
1398 u_c_abrhob*u_c_abnorm_drhoa*c_cab_2 + u_c_ab* &
1399 u_c_abrhobnorm_drhoa*c_cab_2))
1400 e_rb_ndra(ii) = e_rb_ndra(ii) + exc_rhob_norm_drhoa
1401
1402 t1571 = s_anorm_drhoa**2
1403 u_x_anorm_drhoanorm_drhoa = 0.2e1_dp*gamma_x*t1571*t16 - &
1404 0.10e2_dp*t568*t198*t1571 + 0.8e1_dp*t577*t579*t1571
1405 s_a_2norm_drhoanorm_drhoa = 0.2e1_dp*t1571
1406 s_a_21norm_drhoa = s_a_2norm_drhoa
1407 s_avg_2norm_drhoanorm_drhoa = s_a_2norm_drhoanorm_drhoa/0.2e1_dp
1408 s_avg_21norm_drhoa = s_a_21norm_drhoa/0.2e1_dp
1409 t1589 = t336*s_avg_2norm_drhoa
1410 t1590 = t339*s_avg_21norm_drhoa
1411 t1593 = t913*s_avg_2norm_drhoa
1412 u_c_abnorm_drhoanorm_drhoa = gamma_c_ab* &
1413 s_avg_2norm_drhoanorm_drhoa*t162 - 0.2e1_dp*t1589*t1590 + &
1414 0.2e1_dp*t911*t1593*s_avg_21norm_drhoa - t337*t339* &
1415 s_avg_2norm_drhoanorm_drhoa
1416 t1605 = t347*s_a_21norm_drhoa
1417 u_c_anorm_drhoanorm_drhoa = gamma_c_ss*s_a_2norm_drhoanorm_drhoa &
1418 *t165 - 0.2e1_dp*t344*s_a_2norm_drhoa*t1605 + 0.2e1_dp* &
1419 t930*t932*s_a_2norm_drhoa*s_a_21norm_drhoa - t345*t347* &
1420 s_a_2norm_drhoanorm_drhoa
1421
1422 exc_norm_drhoa_norm_drhoa = scale_x*e_lsda_x_a*( &
1423 u_x_anorm_drhoanorm_drhoa*t18 + 0.2e1_dp*u_x_anorm_drhoa**2* &
1424 c_x_2 + u_x_a*u_x_anorm_drhoanorm_drhoa*c_x_2) + scale_c*( &
1425 e_lsda_c_ab*(u_c_abnorm_drhoanorm_drhoa*t170 + 0.2e1_dp* &
1426 u_c_abnorm_drhoa*(gamma_c_ab*s_avg_21norm_drhoa*t162 - t337* &
1427 t1590)*c_cab_2 + u_c_ab*u_c_abnorm_drhoanorm_drhoa*c_cab_2) + &
1428 e_lsda_c_a*(u_c_anorm_drhoanorm_drhoa*t173 + 0.2e1_dp* &
1429 u_c_anorm_drhoa*(gamma_c_ss*s_a_21norm_drhoa*t165 - t345* &
1430 t1605)*c_css_2 + u_c_a*u_c_anorm_drhoanorm_drhoa*c_css_2))
1431 e_ndra_ndra(ii) = e_ndra_ndra(ii) + exc_norm_drhoa_norm_drhoa
1432
1433 u_c_abrhoanorm_drhob = -0.2e1_dp*t906*t539 + 0.2e1_dp*t911* &
1434 t914*s_avg_2norm_drhob
1435
1436 exc_rhoa_norm_drhob = scale_c*(e_lsda_c_abrhoa*gc_abnorm_drhob &
1437 + e_lsda_c_ab*(u_c_abrhoanorm_drhob*t170 + 0.2e1_dp* &
1438 u_c_abrhoa*u_c_abnorm_drhob*c_cab_2 + u_c_ab* &
1439 u_c_abrhoanorm_drhob*c_cab_2))
1440 e_ra_ndrb(ii) = e_ra_ndrb(ii) + exc_rhoa_norm_drhob
1441
1442 s_brhobnorm_drhob = -0.4e1_dp/0.3e1_dp*t367
1443 u_x_brhobnorm_drhob = 0.2e1_dp*gamma_x*s_bnorm_drhob*t371 - &
1444 0.10e2_dp*t1165*t377*s_bnorm_drhob + 0.2e1_dp*t370*t29* &
1445 s_brhobnorm_drhob + 0.8e1_dp*t1173*t1175*s_brhob* &
1446 s_bnorm_drhob - 0.2e1_dp*t374*t376*s_brhobnorm_drhob
1447 s_b_2rhobnorm_drhob = 0.2e1_dp*s_bnorm_drhob*s_brhob + &
1448 0.2e1_dp*s_b*s_brhobnorm_drhob
1449 s_avg_2rhobnorm_drhob = s_b_2rhobnorm_drhob/0.2e1_dp
1450 u_c_abrhobnorm_drhob = gamma_c_ab*s_avg_2rhobnorm_drhob*t162 - &
1451 0.2e1_dp*t1436*t539 + 0.2e1_dp*t911*t1440* &
1452 s_avg_2norm_drhob - t337*t339*s_avg_2rhobnorm_drhob
1453 u_c_brhobnorm_drhob = gamma_c_ss*s_b_2rhobnorm_drhob*t168 - &
1454 0.2e1_dp*t1451*t543 + 0.2e1_dp*t1455*t1458*s_b_2norm_drhob &
1455 - t484*t486*s_b_2rhobnorm_drhob
1456
1457 exc_rhob_norm_drhob = scale_x*(e_lsda_x_brhob*gx_bnorm_drhob + &
1458 e_lsda_x_b*(u_x_brhobnorm_drhob*t31 + 0.2e1_dp*u_x_brhob* &
1459 u_x_bnorm_drhob*c_x_2 + u_x_b*u_x_brhobnorm_drhob*c_x_2)) + &
1460 scale_c*(e_lsda_c_abrhob*gc_abnorm_drhob + e_lsda_c_ab*( &
1461 u_c_abrhobnorm_drhob*t170 + 0.2e1_dp*u_c_abrhob* &
1462 u_c_abnorm_drhob*c_cab_2 + u_c_ab*u_c_abrhobnorm_drhob*c_cab_2 &
1463 ) + e_lsda_c_brhob*gc_bnorm_drhob + e_lsda_c_b*( &
1464 u_c_brhobnorm_drhob*t176 + 0.2e1_dp*u_c_brhob*u_c_bnorm_drhob &
1465 *c_css_2 + u_c_b*u_c_brhobnorm_drhob*c_css_2))
1466 e_rb_ndrb(ii) = e_rb_ndrb(ii) + exc_rhob_norm_drhob
1467
1468 u_c_abnorm_drhoanorm_drhob = -0.2e1_dp*t1589*t539 + 0.2e1_dp* &
1469 t911*t1593*s_avg_2norm_drhob
1470
1471 exc_norm_drhoa_norm_drhob = scale_c*e_lsda_c_ab*( &
1472 u_c_abnorm_drhoanorm_drhob*t170 + 0.2e1_dp*u_c_abnorm_drhoa* &
1473 u_c_abnorm_drhob*c_cab_2 + u_c_ab*u_c_abnorm_drhoanorm_drhob* &
1474 c_cab_2)
1475 e_ndra_ndrb(ii) = e_ndra_ndrb(ii) + exc_norm_drhoa_norm_drhob
1476
1477 t1719 = s_bnorm_drhob**2
1478 u_x_bnorm_drhobnorm_drhob = 0.2e1_dp*gamma_x*t1719*t29 - &
1479 0.10e2_dp*t1165*t376*t1719 + 0.8e1_dp*t1173*t1175*t1719
1480 s_b_2norm_drhobnorm_drhob = 0.2e1_dp*t1719
1481 s_b_21norm_drhob = s_b_2norm_drhob
1482 s_avg_2norm_drhobnorm_drhob = s_b_2norm_drhobnorm_drhob/0.2e1_dp
1483 s_avg_21norm_drhob = s_b_21norm_drhob/0.2e1_dp
1484 t1738 = t339*s_avg_21norm_drhob
1485 u_c_abnorm_drhobnorm_drhob = gamma_c_ab* &
1486 s_avg_2norm_drhobnorm_drhob*t162 - 0.2e1_dp*t336* &
1487 s_avg_2norm_drhob*t1738 + 0.2e1_dp*t911*t913* &
1488 s_avg_2norm_drhob*s_avg_21norm_drhob - t337*t339* &
1489 s_avg_2norm_drhobnorm_drhob
1490 t1753 = t486*s_b_21norm_drhob
1491 u_c_bnorm_drhobnorm_drhob = gamma_c_ss*s_b_2norm_drhobnorm_drhob &
1492 *t168 - 0.2e1_dp*t344*s_b_2norm_drhob*t1753 + 0.2e1_dp* &
1493 t1455*t1457*s_b_2norm_drhob*s_b_21norm_drhob - t484*t486* &
1494 s_b_2norm_drhobnorm_drhob
1495
1496 exc_norm_drhob_norm_drhob = scale_x*e_lsda_x_b*( &
1497 u_x_bnorm_drhobnorm_drhob*t31 + 0.2e1_dp*u_x_bnorm_drhob**2* &
1498 c_x_2 + u_x_b*u_x_bnorm_drhobnorm_drhob*c_x_2) + scale_c*( &
1499 e_lsda_c_ab*(u_c_abnorm_drhobnorm_drhob*t170 + 0.2e1_dp* &
1500 u_c_abnorm_drhob*(gamma_c_ab*s_avg_21norm_drhob*t162 - t337* &
1501 t1738)*c_cab_2 + u_c_ab*u_c_abnorm_drhobnorm_drhob*c_cab_2) + &
1502 e_lsda_c_b*(u_c_bnorm_drhobnorm_drhob*t176 + 0.2e1_dp* &
1503 u_c_bnorm_drhob*(gamma_c_ss*s_b_21norm_drhob*t168 - t484* &
1504 t1753)*c_css_2 + u_c_b*u_c_bnorm_drhobnorm_drhob*c_css_2))
1505 e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii) + exc_norm_drhob_norm_drhob
1506 END IF ! <1 || >1
1507 END IF ! /=0
1508 END IF ! rho>epsilon_rho
1509 END DO
1510!$OMP END DO
1511 END SELECT
1512 END SUBROUTINE b97_lsd_calc
1513
1514! **************************************************************************************************
1515!> \brief low level calculation of the b97 exchange-correlation functional for lda
1516!> \param rho_tot ...
1517!> \param norm_drho || grad (rhoa+rhob) ||
1518!> \param e_0 adds to it the local value of the functional
1519!> \param e_r derivative of the energy density with respect to r
1520!> \param e_r_r derivative of the energy density with respect to r_r
1521!> \param e_ndr derivative of the energy density with respect to ndr
1522!> \param e_r_ndr derivative of the energy density with respect to r_ndr
1523!> \param e_ndr_ndr derivative of the energy density with respect to ndr_ndr
1524!> \param grad_deriv ...
1525!> \param npoints ...
1526!> \param epsilon_rho ...
1527!> \param param ...
1528!> \param scale_c_in ...
1529!> \param scale_x_in ...
1530!> \author fawzi
1531!> \note slow version
1532! **************************************************************************************************
1533 SUBROUTINE b97_lda_calc(rho_tot, norm_drho, &
1534 e_0, e_r, e_r_r, e_ndr, e_r_ndr, e_ndr_ndr, &
1535 grad_deriv, npoints, epsilon_rho, &
1536 param, scale_c_in, scale_x_in)
1537 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho_tot, norm_drho
1538 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_r, e_r_r, e_ndr, e_r_ndr, &
1539 e_ndr_ndr
1540 INTEGER, INTENT(in) :: grad_deriv, npoints
1541 REAL(kind=dp), INTENT(in) :: epsilon_rho
1542 INTEGER, INTENT(in) :: param
1543 REAL(kind=dp), INTENT(in) :: scale_c_in, scale_x_in
1544
1545 INTEGER :: ii
1546 REAL(kind=dp) :: a_1, a_2, a_3, alpha_1_1, alpha_1_2, alpha_1_3, alpha_c, alpha_c1rhoa, &
1547 alpha_c1rhob, alpha_crhoa, alpha_crhob, beta_1_1, beta_1_2, beta_1_3, beta_2_1, beta_2_2, &
1548 beta_2_3, beta_3_1, beta_3_2, beta_3_3, beta_4_1, beta_4_2, beta_4_3, c_cab_0, c_cab_1, &
1549 c_cab_2, c_css_0, c_css_1, c_css_2, c_x_0, c_x_1, c_x_2, chi, chirhoa, chirhoarhoa, &
1550 chirhoarhob, chirhob, chirhobrhob, e_c_u_0, e_c_u_01rhoa, e_c_u_01rhob, e_c_u_0rhoa, &
1551 e_c_u_0rhoarhoa, e_c_u_0rhoarhob, e_c_u_0rhob, e_c_u_0rhobrhob, e_c_u_1rhoa, e_c_u_1rhob, &
1552 e_lsda_c_a, e_lsda_c_a1rhoa, e_lsda_c_ab, e_lsda_c_abrhoa
1553 REAL(kind=dp) :: e_lsda_c_abrhob, e_lsda_c_arhoa, e_lsda_c_arhoarhoa, e_lsda_c_b, &
1554 e_lsda_c_b1rhob, e_lsda_c_brhob, e_lsda_c_brhobrhob, e_lsda_x_a, e_lsda_x_arhoa, &
1555 e_lsda_x_b, e_lsda_x_brhob, epsilon_c_unif, epsilon_c_unif1rhoa, epsilon_c_unif1rhob, &
1556 epsilon_c_unif_a, epsilon_c_unif_a1rhoa, epsilon_c_unif_arhoa, epsilon_c_unif_b, &
1557 epsilon_c_unif_b1rhob, epsilon_c_unif_brhob, epsilon_c_unifrhoa, epsilon_c_unifrhob, exc, &
1558 exc_norm_drhoa, exc_norm_drhoa_norm_drhoa, exc_norm_drhoa_norm_drhob, exc_norm_drhob, &
1559 exc_norm_drhob_norm_drhob, exc_rhoa, exc_rhoa_norm_drhoa, exc_rhoa_norm_drhob
1560 REAL(kind=dp) :: exc_rhoa_rhoa, exc_rhoa_rhob, exc_rhob, exc_rhob_norm_drhoa, &
1561 exc_rhob_norm_drhob, exc_rhob_rhob, f, f1rhoa, f1rhob, frhoa, frhoarhoa, frhoarhob, &
1562 frhob, frhobrhob, gamma_c_ab, gamma_c_ss, gamma_x, gc_a, gc_ab, gc_abnorm_drhoa, &
1563 gc_abnorm_drhob, gc_abrhoa, gc_abrhob, gc_anorm_drhoa, gc_arhoa, gc_b, gc_bnorm_drhob, &
1564 gc_brhob, gx_a, gx_anorm_drhoa, gx_arhoa, gx_b, gx_bnorm_drhob, gx_brhob, my_norm_drhoa, &
1565 my_norm_drhob, my_rhoa, my_rhob, p_1, p_2, p_3, rho, rs, rs_a, rs_arhoa, rs_arhoarhoa, &
1566 rs_b, rs_brhob, rs_brhobrhob, rsrhoa, rsrhoarhoa, rsrhoarhob, rsrhob, rsrhobrhob, s_a
1567 REAL(kind=dp) :: s_a_2, s_a_21norm_drhoa, s_a_21rhoa, s_a_2norm_drhoa, &
1568 s_a_2norm_drhoanorm_drhoa, s_a_2rhoa, s_a_2rhoanorm_drhoa, s_a_2rhoarhoa, s_anorm_drhoa, &
1569 s_arhoa, s_arhoanorm_drhoa, s_arhoarhoa, s_avg_2, s_avg_21norm_drhoa, s_avg_21norm_drhob, &
1570 s_avg_21rhoa, s_avg_21rhob, s_avg_2norm_drhoa, s_avg_2norm_drhoanorm_drhoa, &
1571 s_avg_2norm_drhob, s_avg_2norm_drhobnorm_drhob, s_avg_2rhoa, s_avg_2rhoanorm_drhoa, &
1572 s_avg_2rhoarhoa, s_avg_2rhob, s_avg_2rhobnorm_drhob, s_avg_2rhobrhob, s_b, s_b_2, &
1573 s_b_21norm_drhob, s_b_21rhob, s_b_2norm_drhob, s_b_2norm_drhobnorm_drhob, s_b_2rhob, &
1574 s_b_2rhobnorm_drhob
1575 REAL(kind=dp) :: s_b_2rhobrhob, s_bnorm_drhob, s_brhob, s_brhobnorm_drhob, s_brhobrhob, &
1576 scale_c, scale_x, t1, t101, t1012, t1014, t102, t1047, t1049, t105, t106, t107, t108, &
1577 t110, t1107, t112, t113, t1136, t1152, t1157, t116, t1161, t1165, t117, t1172, t1173, &
1578 t1175, t12, t120, t1201, t1205, t122, t1236, t125, t127, t1270, t128, t129, t1299, t13, &
1579 t1321, t1324, t133, t134, t1341, t1348, t1351, t1360, t1368, t138, t1388, t139, t1394, &
1580 t14, t1406, t1407, t1419, t142, t1422, t1436, t1437, t144, t1440, t1451, t1452, t1455, &
1581 t1457, t1458, t147, t149, t15, t150, t151, t155, t156, t1571, t1589, t1590
1582 REAL(kind=dp) :: t1593, t16, t160, t1605, t161, t162, t163, t164, t165, t166, t167, t168, &
1583 t170, t1719, t173, t1738, t1753, t176, t18, t186, t188, t191, t192, t194, t196, t197, &
1584 t198, t199, t2, t20, t207, t208, t209, t21, t210, t212, t216, t220, t221, t222, t223, &
1585 t224, t228, t232, t235, t236, t237, t239, t243, t244, t245, t246, t25, t250, t256, t257, &
1586 t258, t26, t260, t264, t265, t266, t267, t27, t271, t277, t278, t279, t28, t285, t287, &
1587 t289, t29, t290, t291, t293, t294, t295, t297, t299, t3, t301, t302, t304, t31, t312, &
1588 t313, t314, t315, t316, t320, t324, t327, t328, t33, t336, t337, t338
1589 REAL(kind=dp) :: t339, t34, t344, t345, t346, t347, t35, t36, t365, t367, t37, t370, t371, &
1590 t374, t375, t376, t377, t396, t4, t40, t410, t42, t424, t43, t431, t433, t435, t437, &
1591 t438, t439, t441, t443, t445, t446, t448, t456, t457, t458, t459, t46, t460, t464, t468, &
1592 t471, t472, t48, t480, t484, t485, t486, t49, t50, t51, t512, t516, t539, t543, t55, &
1593 t555, t56, t560, t564, t568, t575, t576, t577, t579, t6, t60, t600, t601, t605, t606, &
1594 t608, t619, t62, t621, t626, t627, t631, t632, t633, t639, t644, t646, t647, t659, t66, &
1595 t661, t662, t663, t67, t671, t673, t678, t679, t68, t683, t689, t69
1596 REAL(kind=dp) :: t694, t7, t707, t709, t710, t711, t719, t721, t726, t727, t73, t731, t737, &
1597 t74, t742, t755, t757, t758, t759, t764, t765, t766, t771, t772, t78, t790, t793, t796, &
1598 t8, t80, t811, t818, t821, t830, t838, t84, t85, t858, t86, t864, t87, t876, t877, t889, &
1599 t892, t906, t907, t91, t911, t913, t914, t92, t925, t926, t929, t930, t932, t933, t94, &
1600 t97, t974, t976, t98, t981, t99, t993, u_c_a, u_c_a1rhoa, u_c_ab, u_c_ab1rhoa, &
1601 u_c_ab1rhob, u_c_abnorm_drhoa, u_c_abnorm_drhoanorm_drhoa, u_c_abnorm_drhoanorm_drhob, &
1602 u_c_abnorm_drhob, u_c_abnorm_drhobnorm_drhob, u_c_abrhoa
1603 REAL(kind=dp) :: u_c_abrhoanorm_drhoa, u_c_abrhoanorm_drhob, u_c_abrhoarhoa, u_c_abrhoarhob, &
1604 u_c_abrhob, u_c_abrhobnorm_drhoa, u_c_abrhobnorm_drhob, u_c_abrhobrhob, u_c_anorm_drhoa, &
1605 u_c_anorm_drhoanorm_drhoa, u_c_arhoa, u_c_arhoanorm_drhoa, u_c_arhoarhoa, u_c_b, &
1606 u_c_b1rhob, u_c_bnorm_drhob, u_c_bnorm_drhobnorm_drhob, u_c_brhob, u_c_brhobnorm_drhob, &
1607 u_c_brhobrhob, u_x_a, u_x_a1rhoa, u_x_anorm_drhoa, u_x_anorm_drhoanorm_drhoa, u_x_arhoa, &
1608 u_x_arhoanorm_drhoa, u_x_arhoarhoa, u_x_b, u_x_b1rhob, u_x_bnorm_drhob, &
1609 u_x_bnorm_drhobnorm_drhob, u_x_brhob, u_x_brhobnorm_drhob, u_x_brhobrhob
1610 REAL(kind=dp), DIMENSION(10) :: coeffs
1611
1612 p_1 = 0.10e1_dp
1613 a_1 = 0.31091e-1_dp
1614 alpha_1_1 = 0.21370e0_dp
1615 beta_1_1 = 0.75957e1_dp
1616 beta_2_1 = 0.35876e1_dp
1617 beta_3_1 = 0.16382e1_dp
1618 beta_4_1 = 0.49294e0_dp
1619 p_2 = 0.10e1_dp
1620 a_2 = 0.15545e-1_dp
1621 alpha_1_2 = 0.20548e0_dp
1622 beta_1_2 = 0.141189e2_dp
1623 beta_2_2 = 0.61977e1_dp
1624 beta_3_2 = 0.33662e1_dp
1625 beta_4_2 = 0.62517e0_dp
1626 p_3 = 0.10e1_dp
1627 a_3 = 0.16887e-1_dp
1628 alpha_1_3 = 0.11125e0_dp
1629 beta_1_3 = 0.10357e2_dp
1630 beta_2_3 = 0.36231e1_dp
1631 beta_3_3 = 0.88026e0_dp
1632 beta_4_3 = 0.49671e0_dp
1633
1634 coeffs = b97_coeffs(param)
1635 c_x_0 = coeffs(1)
1636 c_x_1 = coeffs(2)
1637 c_x_2 = coeffs(3)
1638 c_cab_0 = coeffs(4)
1639 c_cab_1 = coeffs(5)
1640 c_cab_2 = coeffs(6)
1641 c_css_0 = coeffs(7)
1642 c_css_1 = coeffs(8)
1643 c_css_2 = coeffs(9)
1644
1645 scale_c = scale_c_in
1646 scale_x = scale_x_in
1647 IF (scale_x == -1.0_dp) scale_x = coeffs(10)
1648
1649 gamma_x = 0.004_dp
1650 gamma_c_ss = 0.2_dp
1651 gamma_c_ab = 0.006_dp
1652
1653 SELECT CASE (grad_deriv)
1654 CASE default
1655 t1 = 3**(0.1e1_dp/0.3e1_dp)
1656 t2 = 4**(0.1e1_dp/0.3e1_dp)
1657 t3 = t2**2
1658 t4 = t1*t3
1659 t6 = (0.1e1_dp/pi)**(0.1e1_dp/0.3e1_dp)
1660!$OMP DO
1661 DO ii = 1, npoints
1662 my_rhoa = 0.5_dp*max(rho_tot(ii), 0.0_dp)
1663 my_rhob = my_rhoa
1664 rho = my_rhoa + my_rhob
1665 IF (rho > epsilon_rho) THEN
1666 my_rhoa = max(my_rhoa, 0.5_dp*epsilon_rho)
1667 my_rhob = max(my_rhob, 0.5_dp*epsilon_rho)
1668 rho = my_rhoa + my_rhob
1669 my_norm_drhoa = 0.5_dp*norm_drho(ii)
1670 my_norm_drhob = 0.5_dp*norm_drho(ii)
1671 t7 = my_rhoa**(0.1e1_dp/0.3e1_dp)
1672 t8 = t7*my_rhoa
1673 e_lsda_x_a = -0.3e1_dp/0.8e1_dp*t4*t6*t8
1674 t12 = 0.1e1_dp/t8
1675 s_a = my_norm_drhoa*t12
1676 t13 = s_a**2
1677 t14 = gamma_x*t13
1678 t15 = 0.1e1_dp + t14
1679 t16 = 0.1e1_dp/t15
1680 u_x_a = t14*t16
1681 t18 = c_x_1 + u_x_a*c_x_2
1682 gx_a = c_x_0 + u_x_a*t18
1683 t20 = my_rhob**(0.1e1_dp/0.3e1_dp)
1684 t21 = t20*my_rhob
1685 e_lsda_x_b = -0.3e1_dp/0.8e1_dp*t4*t6*t21
1686 t25 = 0.1e1_dp/t21
1687 s_b = my_norm_drhob*t25
1688 t26 = s_b**2
1689 t27 = gamma_x*t26
1690 t28 = 0.1e1_dp + t27
1691 t29 = 0.1e1_dp/t28
1692 u_x_b = t27*t29
1693 t31 = c_x_1 + u_x_b*c_x_2
1694 gx_b = c_x_0 + u_x_b*t31
1695 t33 = my_rhoa - my_rhob
1696 t34 = 0.1e1_dp/rho
1697 chi = t33*t34
1698 t35 = 0.1e1_dp/pi
1699 t36 = t35*t34
1700 t37 = t36**(0.1e1_dp/0.3e1_dp)
1701 rs = t4*t37/0.4e1_dp
1702 t40 = 0.1e1_dp + alpha_1_1*rs
1703 t42 = 0.1e1_dp/a_1
1704 t43 = sqrt(rs)
1705 t46 = t43*rs
1706 t48 = p_1 + 0.1e1_dp
1707 t49 = rs**t48
1708 t50 = beta_4_1*t49
1709 t51 = beta_1_1*t43 + beta_2_1*rs + beta_3_1*t46 + t50
1710 t55 = 0.1e1_dp + t42/t51/0.2e1_dp
1711 t56 = log(t55)
1712 e_c_u_0 = -0.2e1_dp*a_1*t40*t56
1713 t60 = 0.1e1_dp + alpha_1_2*rs
1714 t62 = 0.1e1_dp/a_2
1715 t66 = p_2 + 0.1e1_dp
1716 t67 = rs**t66
1717 t68 = beta_4_2*t67
1718 t69 = beta_1_2*t43 + beta_2_2*rs + beta_3_2*t46 + t68
1719 t73 = 0.1e1_dp + t62/t69/0.2e1_dp
1720 t74 = log(t73)
1721 t78 = 0.1e1_dp + alpha_1_3*rs
1722 t80 = 0.1e1_dp/a_3
1723 t84 = p_3 + 0.1e1_dp
1724 t85 = rs**t84
1725 t86 = beta_4_3*t85
1726 t87 = beta_1_3*t43 + beta_2_3*rs + beta_3_3*t46 + t86
1727 t91 = 0.1e1_dp + t80/t87/0.2e1_dp
1728 t92 = log(t91)
1729 alpha_c = 0.2e1_dp*a_3*t78*t92
1730 t94 = 2**(0.1e1_dp/0.3e1_dp)
1731 t97 = 1/(2*t94 - 2)
1732 t98 = 0.1e1_dp + chi
1733 t99 = t98**(0.1e1_dp/0.3e1_dp)
1734 t101 = 0.1e1_dp - chi
1735 t102 = t101**(0.1e1_dp/0.3e1_dp)
1736 f = (t99*t98 + t102*t101 - 0.2e1_dp)*t97
1737 t105 = alpha_c*f
1738 t106 = 0.9e1_dp/0.8e1_dp/t97
1739 t107 = chi**2
1740 t108 = t107**2
1741 t110 = t106*(0.1e1_dp - t108)
1742 t112 = -0.2e1_dp*a_2*t60*t74 - e_c_u_0
1743 t113 = t112*f
1744 epsilon_c_unif = e_c_u_0 + t105*t110 + t113*t108
1745 t116 = t35/my_rhoa
1746 t117 = t116**(0.1e1_dp/0.3e1_dp)
1747 rs_a = t4*t117/0.4e1_dp
1748 t120 = 0.1e1_dp + alpha_1_2*rs_a
1749 t122 = sqrt(rs_a)
1750 t125 = t122*rs_a
1751 t127 = rs_a**t66
1752 t128 = beta_4_2*t127
1753 t129 = beta_1_2*t122 + beta_2_2*rs_a + beta_3_2*t125 + t128
1754 t133 = 0.1e1_dp + t62/t129/0.2e1_dp
1755 t134 = log(t133)
1756 epsilon_c_unif_a = -0.2e1_dp*a_2*t120*t134
1757 t138 = t35/my_rhob
1758 t139 = t138**(0.1e1_dp/0.3e1_dp)
1759 rs_b = t4*t139/0.4e1_dp
1760 t142 = 0.1e1_dp + alpha_1_2*rs_b
1761 t144 = sqrt(rs_b)
1762 t147 = t144*rs_b
1763 t149 = rs_b**t66
1764 t150 = beta_4_2*t149
1765 t151 = beta_1_2*t144 + beta_2_2*rs_b + beta_3_2*t147 + t150
1766 t155 = 0.1e1_dp + t62/t151/0.2e1_dp
1767 t156 = log(t155)
1768 epsilon_c_unif_b = -0.2e1_dp*a_2*t142*t156
1769 s_a_2 = t13
1770 s_b_2 = t26
1771 s_avg_2 = s_a_2/0.2e1_dp + s_b_2/0.2e1_dp
1772 e_lsda_c_a = epsilon_c_unif_a*my_rhoa
1773 e_lsda_c_b = epsilon_c_unif_b*my_rhob
1774 t160 = gamma_c_ab*s_avg_2
1775 t161 = 0.1e1_dp + t160
1776 t162 = 0.1e1_dp/t161
1777 u_c_ab = t160*t162
1778 t163 = gamma_c_ss*s_a_2
1779 t164 = 0.1e1_dp + t163
1780 t165 = 0.1e1_dp/t164
1781 u_c_a = t163*t165
1782 t166 = gamma_c_ss*s_b_2
1783 t167 = 0.1e1_dp + t166
1784 t168 = 0.1e1_dp/t167
1785 u_c_b = t166*t168
1786 e_lsda_c_ab = epsilon_c_unif*rho - e_lsda_c_a - e_lsda_c_b
1787 t170 = c_cab_1 + u_c_ab*c_cab_2
1788 gc_ab = c_cab_0 + u_c_ab*t170
1789 t173 = c_css_1 + u_c_a*c_css_2
1790 gc_a = c_css_0 + u_c_a*t173
1791 t176 = c_css_1 + u_c_b*c_css_2
1792 gc_b = c_css_0 + u_c_b*t176
1793
1794 exc = scale_x*(e_lsda_x_a*gx_a + e_lsda_x_b*gx_b) + scale_c &
1795 *(e_lsda_c_ab*gc_ab + e_lsda_c_a*gc_a + e_lsda_c_b*gc_b)
1796 e_0(ii) = e_0(ii) + exc
1797
1798 IF (grad_deriv /= 0) THEN
1799
1800 e_lsda_x_arhoa = -t4*t6*t7/0.2e1_dp
1801 t186 = my_rhoa**2
1802 t188 = 0.1e1_dp/t7/t186
1803 s_arhoa = -0.4e1_dp/0.3e1_dp*my_norm_drhoa*t188
1804 t191 = gamma_x*s_a
1805 t192 = t16*s_arhoa
1806 t194 = gamma_x**2
1807 t196 = t194*s_a_2*s_a
1808 t197 = t15**2
1809 t198 = 0.1e1_dp/t197
1810 t199 = t198*s_arhoa
1811 u_x_arhoa = 0.2e1_dp*t191*t192 - 0.2e1_dp*t196*t199
1812 gx_arhoa = u_x_arhoa*t18 + u_x_a*u_x_arhoa*c_x_2
1813 t207 = rho**2
1814 t208 = 0.1e1_dp/t207
1815 t209 = t33*t208
1816 chirhoa = t34 - t209
1817 t210 = t37**2
1818 t212 = 0.1e1_dp/t210*t35
1819 rsrhoa = -t4*t212*t208/0.12e2_dp
1820 t216 = a_1*alpha_1_1
1821 t220 = t51**2
1822 t221 = 0.1e1_dp/t220
1823 t222 = t40*t221
1824 t223 = 0.1e1_dp/t43
1825 t224 = beta_1_1*t223
1826 t228 = beta_3_1*t43
1827 t232 = 0.1e1_dp/rs
1828 t235 = t224*rsrhoa/0.2e1_dp + beta_2_1*rsrhoa + 0.3e1_dp/ &
1829 0.2e1_dp*t228*rsrhoa + t50*t48*rsrhoa*t232
1830 t236 = 0.1e1_dp/t55
1831 t237 = t235*t236
1832 e_c_u_0rhoa = -0.2e1_dp*t216*rsrhoa*t56 + t222*t237
1833 t239 = a_2*alpha_1_2
1834 t243 = t69**2
1835 t244 = 0.1e1_dp/t243
1836 t245 = t60*t244
1837 t246 = beta_1_2*t223
1838 t250 = beta_3_2*t43
1839 t256 = t246*rsrhoa/0.2e1_dp + beta_2_2*rsrhoa + 0.3e1_dp/ &
1840 0.2e1_dp*t250*rsrhoa + t68*t66*rsrhoa*t232
1841 t257 = 0.1e1_dp/t73
1842 t258 = t256*t257
1843 e_c_u_1rhoa = -0.2e1_dp*t239*rsrhoa*t74 + t245*t258
1844 t260 = a_3*alpha_1_3
1845 t264 = t87**2
1846 t265 = 0.1e1_dp/t264
1847 t266 = t78*t265
1848 t267 = beta_1_3*t223
1849 t271 = beta_3_3*t43
1850 t277 = t267*rsrhoa/0.2e1_dp + beta_2_3*rsrhoa + 0.3e1_dp/ &
1851 0.2e1_dp*t271*rsrhoa + t86*t84*rsrhoa*t232
1852 t278 = 0.1e1_dp/t91
1853 t279 = t277*t278
1854 alpha_crhoa = 0.2e1_dp*t260*rsrhoa*t92 - t266*t279
1855 frhoa = (0.4e1_dp/0.3e1_dp*t99*chirhoa - 0.4e1_dp/0.3e1_dp &
1856 *t102*chirhoa)*t97
1857 t285 = alpha_crhoa*f
1858 t287 = alpha_c*frhoa
1859 t289 = t107*chi
1860 t290 = t106*t289
1861 t291 = t290*chirhoa
1862 t293 = 0.4e1_dp*t105*t291
1863 t294 = e_c_u_1rhoa - e_c_u_0rhoa
1864 t295 = t294*f
1865 t297 = t112*frhoa
1866 t299 = t289*chirhoa
1867 t301 = 0.4e1_dp*t113*t299
1868 epsilon_c_unifrhoa = e_c_u_0rhoa + t285*t110 + t287*t110 - &
1869 t293 + t295*t108 + t297*t108 + t301
1870 t302 = t117**2
1871 t304 = 0.1e1_dp/t302*t35
1872 rs_arhoa = -t4*t304/t186/0.12e2_dp
1873 t312 = t129**2
1874 t313 = 0.1e1_dp/t312
1875 t314 = t120*t313
1876 t315 = 0.1e1_dp/t122
1877 t316 = beta_1_2*t315
1878 t320 = beta_3_2*t122
1879 t324 = 0.1e1_dp/rs_a
1880 t327 = t316*rs_arhoa/0.2e1_dp + beta_2_2*rs_arhoa + 0.3e1_dp &
1881 /0.2e1_dp*t320*rs_arhoa + t128*t66*rs_arhoa*t324
1882 t328 = 0.1e1_dp/t133
1883 epsilon_c_unif_arhoa = -0.2e1_dp*t239*rs_arhoa*t134 + t314* &
1884 t327*t328
1885 s_a_2rhoa = 0.2e1_dp*s_a*s_arhoa
1886 s_avg_2rhoa = s_a_2rhoa/0.2e1_dp
1887 e_lsda_c_arhoa = epsilon_c_unif_arhoa*my_rhoa + epsilon_c_unif_a
1888 t336 = gamma_c_ab**2
1889 t337 = t336*s_avg_2
1890 t338 = t161**2
1891 t339 = 0.1e1_dp/t338
1892 u_c_abrhoa = gamma_c_ab*s_avg_2rhoa*t162 - t337*t339*s_avg_2rhoa
1893 t344 = gamma_c_ss**2
1894 t345 = t344*s_a_2
1895 t346 = t164**2
1896 t347 = 0.1e1_dp/t346
1897 u_c_arhoa = gamma_c_ss*s_a_2rhoa*t165 - t345*t347*s_a_2rhoa
1898 e_lsda_c_abrhoa = epsilon_c_unifrhoa*rho + epsilon_c_unif - &
1899 e_lsda_c_arhoa
1900 gc_abrhoa = u_c_abrhoa*t170 + u_c_ab*u_c_abrhoa*c_cab_2
1901 gc_arhoa = u_c_arhoa*t173 + u_c_a*u_c_arhoa*c_css_2
1902
1903 exc_rhoa = scale_x*(e_lsda_x_arhoa*gx_a + e_lsda_x_a* &
1904 gx_arhoa) + scale_c*(e_lsda_c_abrhoa*gc_ab + e_lsda_c_ab* &
1905 gc_abrhoa + e_lsda_c_arhoa*gc_a + e_lsda_c_a*gc_arhoa)
1906 e_r(ii) = e_r(ii) + 0.5_dp*exc_rhoa
1907
1908 e_lsda_x_brhob = -t4*t6*t20/0.2e1_dp
1909 t365 = my_rhob**2
1910 t367 = 0.1e1_dp/t20/t365
1911 s_brhob = -0.4e1_dp/0.3e1_dp*my_norm_drhob*t367
1912 t370 = gamma_x*s_b
1913 t371 = t29*s_brhob
1914 t374 = t194*s_b_2*s_b
1915 t375 = t28**2
1916 t376 = 0.1e1_dp/t375
1917 t377 = t376*s_brhob
1918 u_x_brhob = 0.2e1_dp*t370*t371 - 0.2e1_dp*t374*t377
1919 gx_brhob = u_x_brhob*t31 + u_x_b*u_x_brhob*c_x_2
1920 chirhob = -t34 - t209
1921 rsrhob = rsrhoa
1922 t396 = t224*rsrhob/0.2e1_dp + beta_2_1*rsrhob + 0.3e1_dp/ &
1923 0.2e1_dp*t228*rsrhob + t50*t48*rsrhob*t232
1924 e_c_u_0rhob = -0.2e1_dp*t216*rsrhob*t56 + t222*t396*t236
1925 t410 = t246*rsrhob/0.2e1_dp + beta_2_2*rsrhob + 0.3e1_dp/ &
1926 0.2e1_dp*t250*rsrhob + t68*t66*rsrhob*t232
1927 e_c_u_1rhob = -0.2e1_dp*t239*rsrhob*t74 + t245*t410*t257
1928 t424 = t267*rsrhob/0.2e1_dp + beta_2_3*rsrhob + 0.3e1_dp/ &
1929 0.2e1_dp*t271*rsrhob + t86*t84*rsrhob*t232
1930 alpha_crhob = 0.2e1_dp*t260*rsrhob*t92 - t266*t424*t278
1931 frhob = (0.4e1_dp/0.3e1_dp*t99*chirhob - 0.4e1_dp/0.3e1_dp &
1932 *t102*chirhob)*t97
1933 t431 = alpha_crhob*f
1934 t433 = alpha_c*frhob
1935 t435 = t290*chirhob
1936 t437 = 0.4e1_dp*t105*t435
1937 t438 = e_c_u_1rhob - e_c_u_0rhob
1938 t439 = t438*f
1939 t441 = t112*frhob
1940 t443 = t289*chirhob
1941 t445 = 0.4e1_dp*t113*t443
1942 epsilon_c_unifrhob = e_c_u_0rhob + t431*t110 + t433*t110 - &
1943 t437 + t439*t108 + t441*t108 + t445
1944 t446 = t139**2
1945 t448 = 0.1e1_dp/t446*t35
1946 rs_brhob = -t4*t448/t365/0.12e2_dp
1947 t456 = t151**2
1948 t457 = 0.1e1_dp/t456
1949 t458 = t142*t457
1950 t459 = 0.1e1_dp/t144
1951 t460 = beta_1_2*t459
1952 t464 = beta_3_2*t144
1953 t468 = 0.1e1_dp/rs_b
1954 t471 = t460*rs_brhob/0.2e1_dp + beta_2_2*rs_brhob + 0.3e1_dp &
1955 /0.2e1_dp*rs_brhob*t464 + t150*t66*rs_brhob*t468
1956 t472 = 0.1e1_dp/t155
1957 epsilon_c_unif_brhob = -0.2e1_dp*t239*rs_brhob*t156 + t458* &
1958 t471*t472
1959 s_b_2rhob = 0.2e1_dp*s_b*s_brhob
1960 s_avg_2rhob = s_b_2rhob/0.2e1_dp
1961 e_lsda_c_brhob = epsilon_c_unif_brhob*my_rhob + epsilon_c_unif_b
1962 t480 = t339*s_avg_2rhob
1963 u_c_abrhob = gamma_c_ab*s_avg_2rhob*t162 - t337*t480
1964 t484 = t344*s_b_2
1965 t485 = t167**2
1966 t486 = 0.1e1_dp/t485
1967 u_c_brhob = gamma_c_ss*s_b_2rhob*t168 - t484*t486*s_b_2rhob
1968 e_lsda_c_abrhob = epsilon_c_unifrhob*rho + epsilon_c_unif - &
1969 e_lsda_c_brhob
1970 gc_abrhob = u_c_abrhob*t170 + u_c_ab*u_c_abrhob*c_cab_2
1971 gc_brhob = u_c_brhob*t176 + u_c_b*u_c_brhob*c_css_2
1972
1973 exc_rhob = scale_x*(e_lsda_x_brhob*gx_b + e_lsda_x_b* &
1974 gx_brhob) + scale_c*(e_lsda_c_abrhob*gc_ab + e_lsda_c_ab* &
1975 gc_abrhob + e_lsda_c_brhob*gc_b + e_lsda_c_b*gc_brhob)
1976 e_r(ii) = e_r(ii) + 0.5_dp*exc_rhob
1977
1978 s_anorm_drhoa = t12
1979 u_x_anorm_drhoa = 0.2e1_dp*t191*t16*s_anorm_drhoa - 0.2e1_dp &
1980 *t196*t198*s_anorm_drhoa
1981 gx_anorm_drhoa = u_x_anorm_drhoa*t18 + u_x_a*u_x_anorm_drhoa*c_x_2
1982 s_a_2norm_drhoa = 0.2e1_dp*s_a*s_anorm_drhoa
1983 s_avg_2norm_drhoa = s_a_2norm_drhoa/0.2e1_dp
1984 t512 = t339*s_avg_2norm_drhoa
1985 u_c_abnorm_drhoa = gamma_c_ab*s_avg_2norm_drhoa*t162 - t337*t512
1986 t516 = t347*s_a_2norm_drhoa
1987 u_c_anorm_drhoa = gamma_c_ss*s_a_2norm_drhoa*t165 - t345*t516
1988 gc_abnorm_drhoa = u_c_abnorm_drhoa*t170 + u_c_ab* &
1989 u_c_abnorm_drhoa*c_cab_2
1990 gc_anorm_drhoa = u_c_anorm_drhoa*t173 + u_c_a*u_c_anorm_drhoa &
1991 *c_css_2
1992
1993 exc_norm_drhoa = scale_x*e_lsda_x_a*gx_anorm_drhoa + scale_c* &
1994 (e_lsda_c_ab*gc_abnorm_drhoa + e_lsda_c_a*gc_anorm_drhoa)
1995 e_ndr(ii) = e_ndr(ii) + 0.5_dp*exc_norm_drhoa
1996
1997 s_bnorm_drhob = t25
1998 u_x_bnorm_drhob = 0.2e1_dp*t370*t29*s_bnorm_drhob - 0.2e1_dp &
1999 *t374*t376*s_bnorm_drhob
2000 gx_bnorm_drhob = u_x_bnorm_drhob*t31 + u_x_b*u_x_bnorm_drhob*c_x_2
2001 s_b_2norm_drhob = 0.2e1_dp*s_b*s_bnorm_drhob
2002 s_avg_2norm_drhob = s_b_2norm_drhob/0.2e1_dp
2003 t539 = t339*s_avg_2norm_drhob
2004 u_c_abnorm_drhob = gamma_c_ab*s_avg_2norm_drhob*t162 - t337*t539
2005 t543 = t486*s_b_2norm_drhob
2006 u_c_bnorm_drhob = gamma_c_ss*s_b_2norm_drhob*t168 - t484*t543
2007 gc_abnorm_drhob = u_c_abnorm_drhob*t170 + u_c_ab* &
2008 u_c_abnorm_drhob*c_cab_2
2009 gc_bnorm_drhob = u_c_bnorm_drhob*t176 + u_c_b*u_c_bnorm_drhob &
2010 *c_css_2
2011
2012 exc_norm_drhob = scale_x*e_lsda_x_b*gx_bnorm_drhob + scale_c* &
2013 (e_lsda_c_ab*gc_abnorm_drhob + e_lsda_c_b*gc_bnorm_drhob)
2014 e_ndr(ii) = e_ndr(ii) + 0.5_dp*exc_norm_drhob
2015
2016 IF (grad_deriv > 1) THEN
2017 t555 = t7**2
2018 t560 = t186*my_rhoa
2019 s_arhoarhoa = 0.28e2_dp/0.9e1_dp*my_norm_drhoa/t7/t560
2020 t564 = s_arhoa**2
2021 t568 = t194*s_a_2
2022 t575 = t194*gamma_x
2023 t576 = s_a_2**2
2024 t577 = t575*t576
2025 t579 = 0.1e1_dp/t197/t15
2026 u_x_arhoarhoa = 0.2e1_dp*gamma_x*t564*t16 - 0.10e2_dp*t568 &
2027 *t198*t564 + 0.2e1_dp*t191*t16*s_arhoarhoa + 0.8e1_dp* &
2028 t577*t579*t564 - 0.2e1_dp*t196*t198*s_arhoarhoa
2029 u_x_a1rhoa = u_x_arhoa
2030 t600 = 0.1e1_dp/t207/rho
2031 t601 = t33*t600
2032 chirhoarhoa = -0.2e1_dp*t208 + 0.2e1_dp*t601
2033 t605 = 0.3141592654e1_dp**2
2034 t606 = 0.1e1_dp/t605
2035 t608 = t207**2
2036 rsrhoarhoa = -t4/t210/t36*t606/t608/0.18e2_dp + &
2037 t4*t212*t600/0.6e1_dp
2038 t619 = alpha_1_1*rsrhoa
2039 t621 = t221*t235*t236
2040 t626 = t40/t220/t51
2041 t627 = t235**2
2042 t631 = 0.1e1_dp/t46
2043 t632 = beta_1_1*t631
2044 t633 = rsrhoa**2
2045 t639 = beta_3_1*t223
2046 t644 = t48**2
2047 t646 = rs**2
2048 t647 = 0.1e1_dp/t646
2049 t659 = t220**2
2050 t661 = t40/t659
2051 t662 = t55**2
2052 t663 = 0.1e1_dp/t662
2053 e_c_u_0rhoarhoa = -0.2e1_dp*t216*rsrhoarhoa*t56 + 0.2e1_dp* &
2054 t619*t621 - 0.2e1_dp*t626*t627*t236 + t222*(-t632*t633 &
2055 /0.4e1_dp + t224*rsrhoarhoa/0.2e1_dp + beta_2_1*rsrhoarhoa + &
2056 0.3e1_dp/0.4e1_dp*t639*t633 + 0.3e1_dp/0.2e1_dp*t228* &
2057 rsrhoarhoa + t50*t644*t633*t647 + t50*t48*rsrhoarhoa* &
2058 t232 - t50*t48*t633*t647)*t236 + t661*t627*t663*t42/ &
2059 0.2e1_dp
2060 e_c_u_01rhoa = e_c_u_0rhoa
2061 t671 = alpha_1_2*rsrhoa
2062 t673 = t244*t256*t257
2063 t678 = t60/t243/t69
2064 t679 = t256**2
2065 t683 = beta_1_2*t631
2066 t689 = beta_3_2*t223
2067 t694 = t66**2
2068 t707 = t243**2
2069 t709 = t60/t707
2070 t710 = t73**2
2071 t711 = 0.1e1_dp/t710
2072 t719 = alpha_1_3*rsrhoa
2073 t721 = t265*t277*t278
2074 t726 = t78/t264/t87
2075 t727 = t277**2
2076 t731 = beta_1_3*t631
2077 t737 = beta_3_3*t223
2078 t742 = t84**2
2079 t755 = t264**2
2080 t757 = t78/t755
2081 t758 = t91**2
2082 t759 = 0.1e1_dp/t758
2083 alpha_c1rhoa = alpha_crhoa
2084 t764 = t99**2
2085 t765 = 0.1e1_dp/t764
2086 t766 = chirhoa**2
2087 t771 = t102**2
2088 t772 = 0.1e1_dp/t771
2089 frhoarhoa = (0.4e1_dp/0.9e1_dp*t765*t766 + 0.4e1_dp/ &
2090 0.3e1_dp*t99*chirhoarhoa + 0.4e1_dp/0.9e1_dp*t772*t766 - &
2091 0.4e1_dp/0.3e1_dp*t102*chirhoarhoa)*t97
2092 f1rhoa = frhoa
2093 t790 = alpha_c1rhoa*f
2094 t793 = alpha_c*f1rhoa
2095 t796 = t106*t107
2096 t811 = e_c_u_1rhoa - e_c_u_01rhoa
2097 t818 = t811*f
2098 t821 = t112*f1rhoa
2099 t830 = -0.4e1_dp*t105*t290*chirhoarhoa + (-0.2e1_dp*t239* &
2100 rsrhoarhoa*t74 + 0.2e1_dp*t671*t673 - 0.2e1_dp*t678*t679 &
2101 *t257 + t245*(-t683*t633/0.4e1_dp + t246*rsrhoarhoa/ &
2102 0.2e1_dp + beta_2_2*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t689*t633 &
2103 + 0.3e1_dp/0.2e1_dp*t250*rsrhoarhoa + t68*t694*t633* &
2104 t647 + t68*t66*rsrhoarhoa*t232 - t68*t66*t633*t647)* &
2105 t257 + t709*t679*t711*t62/0.2e1_dp - e_c_u_0rhoarhoa)*f* &
2106 t108 + t294*f1rhoa*t108 + 0.4e1_dp*t295*t299 + t811*frhoa &
2107 *t108 + t112*frhoarhoa*t108 + 0.4e1_dp*t297*t299 + 0.4e1_dp &
2108 *t818*t299 + 0.4e1_dp*t821*t299 + 0.12e2_dp*t113*t107* &
2109 t766 + 0.4e1_dp*t113*t289*chirhoarhoa
2110 epsilon_c_unif1rhoa = e_c_u_01rhoa + t790*t110 + t793*t110 - &
2111 t293 + t818*t108 + t821*t108 + t301
2112 t838 = t186**2
2113 rs_arhoarhoa = -t4/t302/t116*t606/t838/0.18e2_dp + &
2114 t4*t304/t560/0.6e1_dp
2115 t858 = t327**2
2116 t864 = rs_arhoa**2
2117 t876 = rs_a**2
2118 t877 = 0.1e1_dp/t876
2119 t889 = t312**2
2120 t892 = t133**2
2121 epsilon_c_unif_a1rhoa = epsilon_c_unif_arhoa
2122 s_a_2rhoarhoa = 0.2e1_dp*t564 + 0.2e1_dp*s_a*s_arhoarhoa
2123 s_a_21rhoa = s_a_2rhoa
2124 s_avg_2rhoarhoa = s_a_2rhoarhoa/0.2e1_dp
2125 s_avg_21rhoa = s_a_21rhoa/0.2e1_dp
2126 e_lsda_c_arhoarhoa = (-0.2e1_dp*t239*rs_arhoarhoa*t134 + &
2127 0.2e1_dp*alpha_1_2*rs_arhoa*t313*t327*t328 - 0.2e1_dp* &
2128 t120/t312/t129*t858*t328 + t314*(-beta_1_2/t125*t864/ &
2129 0.4e1_dp + t316*rs_arhoarhoa/0.2e1_dp + beta_2_2*rs_arhoarhoa &
2130 + 0.3e1_dp/0.4e1_dp*beta_3_2*t315*t864 + 0.3e1_dp/ &
2131 0.2e1_dp*t320*rs_arhoarhoa + t128*t694*t864*t877 + t128* &
2132 t66*rs_arhoarhoa*t324 - t128*t66*t864*t877)*t328 + t120 &
2133 /t889*t858/t892*t62/0.2e1_dp)*my_rhoa + epsilon_c_unif_arhoa &
2134 + epsilon_c_unif_a1rhoa
2135 e_lsda_c_a1rhoa = epsilon_c_unif_a1rhoa*my_rhoa + epsilon_c_unif_a
2136 t906 = t336*s_avg_2rhoa
2137 t907 = t339*s_avg_21rhoa
2138 t911 = t336*gamma_c_ab*s_avg_2
2139 t913 = 0.1e1_dp/t338/t161
2140 t914 = t913*s_avg_2rhoa
2141 u_c_abrhoarhoa = gamma_c_ab*s_avg_2rhoarhoa*t162 - 0.2e1_dp* &
2142 t906*t907 + 0.2e1_dp*t911*t914*s_avg_21rhoa - t337*t339* &
2143 s_avg_2rhoarhoa
2144 u_c_ab1rhoa = gamma_c_ab*s_avg_21rhoa*t162 - t337*t907
2145 t925 = t344*s_a_2rhoa
2146 t926 = t347*s_a_21rhoa
2147 t929 = t344*gamma_c_ss
2148 t930 = t929*s_a_2
2149 t932 = 0.1e1_dp/t346/t164
2150 t933 = t932*s_a_2rhoa
2151 u_c_arhoarhoa = gamma_c_ss*s_a_2rhoarhoa*t165 - 0.2e1_dp* &
2152 t925*t926 + 0.2e1_dp*t930*t933*s_a_21rhoa - t345*t347* &
2153 s_a_2rhoarhoa
2154 u_c_a1rhoa = gamma_c_ss*s_a_21rhoa*t165 - t345*t926
2155
2156 exc_rhoa_rhoa = scale_x*(-t4*t6/t555*gx_a/0.6e1_dp &
2157 + e_lsda_x_arhoa*(u_x_a1rhoa*t18 + u_x_a*u_x_a1rhoa*c_x_2) &
2158 + e_lsda_x_arhoa*gx_arhoa + e_lsda_x_a*(u_x_arhoarhoa*t18 + &
2159 0.2e1_dp*u_x_arhoa*u_x_a1rhoa*c_x_2 + u_x_a*u_x_arhoarhoa* &
2160 c_x_2)) + scale_c*(((e_c_u_0rhoarhoa + (0.2e1_dp*t260* &
2161 rsrhoarhoa*t92 - 0.2e1_dp*t719*t721 + 0.2e1_dp*t726*t727* &
2162 t278 - t266*(-t731*t633/0.4e1_dp + t267*rsrhoarhoa/ &
2163 0.2e1_dp + beta_2_3*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t737*t633 &
2164 + 0.3e1_dp/0.2e1_dp*t271*rsrhoarhoa + t86*t742*t633* &
2165 t647 + t86*t84*rsrhoarhoa*t232 - t86*t84*t633*t647)* &
2166 t278 - t757*t727*t759*t80/0.2e1_dp)*f*t110 + alpha_crhoa &
2167 *f1rhoa*t110 - 0.4e1_dp*t285*t291 + alpha_c1rhoa*frhoa* &
2168 t110 + alpha_c*frhoarhoa*t110 - 0.4e1_dp*t287*t291 - &
2169 0.4e1_dp*t790*t291 - 0.4e1_dp*t793*t291 - 0.12e2_dp*t105* &
2170 t796*t766 + t830)*rho + epsilon_c_unifrhoa + &
2171 epsilon_c_unif1rhoa - e_lsda_c_arhoarhoa)*gc_ab + e_lsda_c_abrhoa &
2172 *(u_c_ab1rhoa*t170 + u_c_ab*u_c_ab1rhoa*c_cab_2) + ( &
2173 epsilon_c_unif1rhoa*rho + epsilon_c_unif - e_lsda_c_a1rhoa)* &
2174 gc_abrhoa + e_lsda_c_ab*(u_c_abrhoarhoa*t170 + 0.2e1_dp* &
2175 u_c_abrhoa*u_c_ab1rhoa*c_cab_2 + u_c_ab*u_c_abrhoarhoa* &
2176 c_cab_2) + e_lsda_c_arhoarhoa*gc_a + e_lsda_c_arhoa*(u_c_a1rhoa &
2177 *t173 + u_c_a*u_c_a1rhoa*c_css_2) + e_lsda_c_a1rhoa*gc_arhoa &
2178 + e_lsda_c_a*(u_c_arhoarhoa*t173 + 0.2e1_dp*u_c_arhoa* &
2179 u_c_a1rhoa*c_css_2 + u_c_a*u_c_arhoarhoa*c_css_2))
2180 e_r_r(ii) = e_r_r(ii) + 0.5_dp*0.5_dp*exc_rhoa_rhoa
2181
2182 chirhoarhob = 0.2e1_dp*t601
2183 rsrhoarhob = rsrhoarhoa
2184 t974 = t221*t396*t236
2185 t976 = alpha_1_1*rsrhob
2186 t981 = rsrhoa*rsrhob
2187 t993 = rsrhob*t647*rsrhoa
2188 e_c_u_0rhoarhob = -0.2e1_dp*t216*rsrhoarhob*t56 + t619* &
2189 t974 + t976*t621 - 0.2e1_dp*t626*t237*t396 + t222*(-t632* &
2190 t981/0.4e1_dp + t224*rsrhoarhob/0.2e1_dp + beta_2_1* &
2191 rsrhoarhob + 0.3e1_dp/0.4e1_dp*t639*t981 + 0.3e1_dp/0.2e1_dp &
2192 *t228*rsrhoarhob + t50*t644*t993 + t50*t48*rsrhoarhob* &
2193 t232 - t50*t48*t993)*t236 + t661*t235*t663*t42*t396/ &
2194 0.2e1_dp
2195 t1012 = t244*t410*t257
2196 t1014 = alpha_1_2*rsrhob
2197 t1047 = t265*t424*t278
2198 t1049 = alpha_1_3*rsrhob
2199 frhoarhob = (0.4e1_dp/0.9e1_dp*t765*chirhoa*chirhob + &
2200 0.4e1_dp/0.3e1_dp*t99*chirhoarhob + 0.4e1_dp/0.9e1_dp*t772 &
2201 *chirhoa*chirhob - 0.4e1_dp/0.3e1_dp*t102*chirhoarhob)* &
2202 t97
2203 t1107 = t107*chirhoa*chirhob
2204 t1136 = -0.4e1_dp*t105*t290*chirhoarhob + (-0.2e1_dp*t239 &
2205 *rsrhoarhob*t74 + t671*t1012 + t1014*t673 - 0.2e1_dp*t678* &
2206 t258*t410 + t245*(-t683*t981/0.4e1_dp + t246*rsrhoarhob/ &
2207 0.2e1_dp + beta_2_2*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t689* &
2208 t981 + 0.3e1_dp/0.2e1_dp*t250*rsrhoarhob + t68*t694*t993 + &
2209 t68*t66*rsrhoarhob*t232 - t68*t66*t993)*t257 + t709* &
2210 t256*t711*t62*t410/0.2e1_dp - e_c_u_0rhoarhob)*f*t108 + &
2211 t294*frhob*t108 + 0.4e1_dp*t295*t443 + t438*frhoa*t108 + &
2212 t112*frhoarhob*t108 + 0.4e1_dp*t297*t443 + 0.4e1_dp*t439 &
2213 *t299 + 0.4e1_dp*t441*t299 + 0.12e2_dp*t113*t1107 + &
2214 0.4e1_dp*t113*t289*chirhoarhob
2215 u_c_abrhoarhob = -0.2e1_dp*t906*t480 + 0.2e1_dp*t911*t914 &
2216 *s_avg_2rhob
2217
2218 exc_rhoa_rhob = scale_c*(((e_c_u_0rhoarhob + (0.2e1_dp*t260* &
2219 rsrhoarhob*t92 - t719*t1047 - t1049*t721 + 0.2e1_dp*t726* &
2220 t279*t424 - t266*(-t731*t981/0.4e1_dp + t267*rsrhoarhob/ &
2221 0.2e1_dp + beta_2_3*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t737*t981 &
2222 + 0.3e1_dp/0.2e1_dp*t271*rsrhoarhob + t86*t742*t993 + t86 &
2223 *t84*rsrhoarhob*t232 - t86*t84*t993)*t278 - t757*t277 &
2224 *t759*t80*t424/0.2e1_dp)*f*t110 + alpha_crhoa*frhob* &
2225 t110 - 0.4e1_dp*t285*t435 + alpha_crhob*frhoa*t110 + alpha_c &
2226 *frhoarhob*t110 - 0.4e1_dp*t287*t435 - 0.4e1_dp*t431* &
2227 t291 - 0.4e1_dp*t433*t291 - 0.12e2_dp*t105*t106*t1107 + &
2228 t1136)*rho + epsilon_c_unifrhoa + epsilon_c_unifrhob)*gc_ab + &
2229 e_lsda_c_abrhoa*gc_abrhob + e_lsda_c_abrhob*gc_abrhoa + &
2230 e_lsda_c_ab*(u_c_abrhoarhob*t170 + 0.2e1_dp*u_c_abrhoa* &
2231 u_c_abrhob*c_cab_2 + u_c_ab*u_c_abrhoarhob*c_cab_2))
2232 e_r_r(ii) = e_r_r(ii) + 0.5_dp*exc_rhoa_rhob
2233
2234 t1152 = t20**2
2235 t1157 = t365*my_rhob
2236 s_brhobrhob = 0.28e2_dp/0.9e1_dp*my_norm_drhob/t20/t1157
2237 t1161 = s_brhob**2
2238 t1165 = t194*s_b_2
2239 t1172 = s_b_2**2
2240 t1173 = t575*t1172
2241 t1175 = 0.1e1_dp/t375/t28
2242 u_x_brhobrhob = 0.2e1_dp*gamma_x*t1161*t29 - 0.10e2_dp* &
2243 t1165*t376*t1161 + 0.2e1_dp*t370*t29*s_brhobrhob + &
2244 0.8e1_dp*t1173*t1175*t1161 - 0.2e1_dp*t374*t376* &
2245 s_brhobrhob
2246 u_x_b1rhob = u_x_brhob
2247 chirhobrhob = 0.2e1_dp*t208 + 0.2e1_dp*t601
2248 rsrhobrhob = rsrhoarhob
2249 t1201 = t396**2
2250 t1205 = rsrhob**2
2251 e_c_u_0rhobrhob = -0.2e1_dp*t216*rsrhobrhob*t56 + 0.2e1_dp* &
2252 t976*t974 - 0.2e1_dp*t626*t1201*t236 + t222*(-t632* &
2253 t1205/0.4e1_dp + t224*rsrhobrhob/0.2e1_dp + beta_2_1* &
2254 rsrhobrhob + 0.3e1_dp/0.4e1_dp*t639*t1205 + 0.3e1_dp/ &
2255 0.2e1_dp*t228*rsrhobrhob + t50*t644*t1205*t647 + t50*t48 &
2256 *rsrhobrhob*t232 - t50*t48*t1205*t647)*t236 + t661* &
2257 t1201*t663*t42/0.2e1_dp
2258 e_c_u_01rhob = e_c_u_0rhob
2259 t1236 = t410**2
2260 t1270 = t424**2
2261 alpha_c1rhob = alpha_crhob
2262 t1299 = chirhob**2
2263 frhobrhob = (0.4e1_dp/0.9e1_dp*t765*t1299 + 0.4e1_dp/ &
2264 0.3e1_dp*t99*chirhobrhob + 0.4e1_dp/0.9e1_dp*t772*t1299 - &
2265 0.4e1_dp/0.3e1_dp*t102*chirhobrhob)*t97
2266 f1rhob = frhob
2267 t1321 = alpha_c1rhob*f
2268 t1324 = alpha_c*f1rhob
2269 t1341 = e_c_u_1rhob - e_c_u_01rhob
2270 t1348 = t1341*f
2271 t1351 = t112*f1rhob
2272 t1360 = -0.4e1_dp*t105*t290*chirhobrhob + (-0.2e1_dp*t239 &
2273 *rsrhobrhob*t74 + 0.2e1_dp*t1014*t1012 - 0.2e1_dp*t678* &
2274 t1236*t257 + t245*(-t683*t1205/0.4e1_dp + t246*rsrhobrhob &
2275 /0.2e1_dp + beta_2_2*rsrhobrhob + 0.3e1_dp/0.4e1_dp*t689* &
2276 t1205 + 0.3e1_dp/0.2e1_dp*t250*rsrhobrhob + t68*t694*t1205 &
2277 *t647 + t68*t66*rsrhobrhob*t232 - t68*t66*t1205*t647) &
2278 *t257 + t709*t1236*t711*t62/0.2e1_dp - e_c_u_0rhobrhob)*f &
2279 *t108 + t438*f1rhob*t108 + 0.4e1_dp*t439*t443 + t1341* &
2280 frhob*t108 + t112*frhobrhob*t108 + 0.4e1_dp*t441*t443 + &
2281 0.4e1_dp*t1348*t443 + 0.4e1_dp*t1351*t443 + 0.12e2_dp*t113 &
2282 *t107*t1299 + 0.4e1_dp*t113*t289*chirhobrhob
2283 epsilon_c_unif1rhob = e_c_u_01rhob + t1321*t110 + t1324*t110 - &
2284 t437 + t1348*t108 + t1351*t108 + t445
2285 t1368 = t365**2
2286 rs_brhobrhob = -t4/t446/t138*t606/t1368/0.18e2_dp &
2287 + t4*t448/t1157/0.6e1_dp
2288 t1388 = t471**2
2289 t1394 = rs_brhob**2
2290 t1406 = rs_b**2
2291 t1407 = 0.1e1_dp/t1406
2292 t1419 = t456**2
2293 t1422 = t155**2
2294 epsilon_c_unif_b1rhob = epsilon_c_unif_brhob
2295 s_b_2rhobrhob = 0.2e1_dp*t1161 + 0.2e1_dp*s_b*s_brhobrhob
2296 s_b_21rhob = s_b_2rhob
2297 s_avg_2rhobrhob = s_b_2rhobrhob/0.2e1_dp
2298 s_avg_21rhob = s_b_21rhob/0.2e1_dp
2299 e_lsda_c_brhobrhob = (-0.2e1_dp*t239*rs_brhobrhob*t156 + &
2300 0.2e1_dp*alpha_1_2*rs_brhob*t457*t471*t472 - 0.2e1_dp* &
2301 t142/t456/t151*t1388*t472 + t458*(-beta_1_2/t147*t1394 &
2302 /0.4e1_dp + t460*rs_brhobrhob/0.2e1_dp + beta_2_2* &
2303 rs_brhobrhob + 0.3e1_dp/0.4e1_dp*beta_3_2*t459*t1394 + &
2304 0.3e1_dp/0.2e1_dp*t464*rs_brhobrhob + t150*t694*t1394* &
2305 t1407 + t150*t66*rs_brhobrhob*t468 - t150*t66*t1394* &
2306 t1407)*t472 + t142/t1419*t1388/t1422*t62/0.2e1_dp)* &
2307 my_rhob + epsilon_c_unif_brhob + epsilon_c_unif_b1rhob
2308 e_lsda_c_b1rhob = epsilon_c_unif_b1rhob*my_rhob + epsilon_c_unif_b
2309 t1436 = t336*s_avg_2rhob
2310 t1437 = t339*s_avg_21rhob
2311 t1440 = t913*s_avg_2rhob
2312 u_c_abrhobrhob = gamma_c_ab*s_avg_2rhobrhob*t162 - 0.2e1_dp* &
2313 t1436*t1437 + 0.2e1_dp*t911*t1440*s_avg_21rhob - t337*t339 &
2314 *s_avg_2rhobrhob
2315 u_c_ab1rhob = gamma_c_ab*s_avg_21rhob*t162 - t337*t1437
2316 t1451 = t344*s_b_2rhob
2317 t1452 = t486*s_b_21rhob
2318 t1455 = t929*s_b_2
2319 t1457 = 0.1e1_dp/t485/t167
2320 t1458 = t1457*s_b_2rhob
2321 u_c_brhobrhob = gamma_c_ss*s_b_2rhobrhob*t168 - 0.2e1_dp* &
2322 t1451*t1452 + 0.2e1_dp*t1455*t1458*s_b_21rhob - t484*t486 &
2323 *s_b_2rhobrhob
2324 u_c_b1rhob = gamma_c_ss*s_b_21rhob*t168 - t484*t1452
2325
2326 exc_rhob_rhob = scale_x*(-t4*t6/t1152*gx_b/ &
2327 0.6e1_dp + e_lsda_x_brhob*(u_x_b1rhob*t31 + u_x_b*u_x_b1rhob* &
2328 c_x_2) + e_lsda_x_brhob*gx_brhob + e_lsda_x_b*(u_x_brhobrhob* &
2329 t31 + 0.2e1_dp*u_x_brhob*u_x_b1rhob*c_x_2 + u_x_b* &
2330 u_x_brhobrhob*c_x_2)) + scale_c*(((e_c_u_0rhobrhob + (0.2e1_dp* &
2331 t260*rsrhobrhob*t92 - 0.2e1_dp*t1049*t1047 + 0.2e1_dp* &
2332 t726*t1270*t278 - t266*(-t731*t1205/0.4e1_dp + t267* &
2333 rsrhobrhob/0.2e1_dp + beta_2_3*rsrhobrhob + 0.3e1_dp/0.4e1_dp &
2334 *t737*t1205 + 0.3e1_dp/0.2e1_dp*t271*rsrhobrhob + t86* &
2335 t742*t1205*t647 + t86*t84*rsrhobrhob*t232 - t86*t84* &
2336 t1205*t647)*t278 - t757*t1270*t759*t80/0.2e1_dp)*f* &
2337 t110 + alpha_crhob*f1rhob*t110 - 0.4e1_dp*t431*t435 + &
2338 alpha_c1rhob*frhob*t110 + alpha_c*frhobrhob*t110 - 0.4e1_dp &
2339 *t433*t435 - 0.4e1_dp*t1321*t435 - 0.4e1_dp*t1324*t435 - &
2340 0.12e2_dp*t105*t796*t1299 + t1360)*rho + epsilon_c_unifrhob &
2341 + epsilon_c_unif1rhob - e_lsda_c_brhobrhob)*gc_ab + &
2342 e_lsda_c_abrhob*(u_c_ab1rhob*t170 + u_c_ab*u_c_ab1rhob* &
2343 c_cab_2) + (epsilon_c_unif1rhob*rho + epsilon_c_unif - &
2344 e_lsda_c_b1rhob)*gc_abrhob + e_lsda_c_ab*(u_c_abrhobrhob*t170 &
2345 + 0.2e1_dp*u_c_abrhob*u_c_ab1rhob*c_cab_2 + u_c_ab* &
2346 u_c_abrhobrhob*c_cab_2) + e_lsda_c_brhobrhob*gc_b + &
2347 e_lsda_c_brhob*(u_c_b1rhob*t176 + u_c_b*u_c_b1rhob*c_css_2) &
2348 + e_lsda_c_b1rhob*gc_brhob + e_lsda_c_b*(u_c_brhobrhob*t176 + &
2349 0.2e1_dp*u_c_brhob*u_c_b1rhob*c_css_2 + u_c_b*u_c_brhobrhob &
2350 *c_css_2))
2351 e_r_r(ii) = e_r_r(ii) + 0.5_dp*0.5_dp*exc_rhob_rhob
2352
2353 s_arhoanorm_drhoa = -0.4e1_dp/0.3e1_dp*t188
2354 u_x_arhoanorm_drhoa = 0.2e1_dp*gamma_x*s_anorm_drhoa*t192 - &
2355 0.10e2_dp*t568*t199*s_anorm_drhoa + 0.2e1_dp*t191*t16* &
2356 s_arhoanorm_drhoa + 0.8e1_dp*t577*t579*s_arhoa*s_anorm_drhoa &
2357 - 0.2e1_dp*t196*t198*s_arhoanorm_drhoa
2358 s_a_2rhoanorm_drhoa = 0.2e1_dp*s_anorm_drhoa*s_arhoa + &
2359 0.2e1_dp*s_a*s_arhoanorm_drhoa
2360 s_avg_2rhoanorm_drhoa = s_a_2rhoanorm_drhoa/0.2e1_dp
2361 u_c_abrhoanorm_drhoa = gamma_c_ab*s_avg_2rhoanorm_drhoa*t162 - &
2362 0.2e1_dp*t906*t512 + 0.2e1_dp*t911*t914*s_avg_2norm_drhoa &
2363 - t337*t339*s_avg_2rhoanorm_drhoa
2364 u_c_arhoanorm_drhoa = gamma_c_ss*s_a_2rhoanorm_drhoa*t165 - &
2365 0.2e1_dp*t925*t516 + 0.2e1_dp*t930*t933*s_a_2norm_drhoa - &
2366 t345*t347*s_a_2rhoanorm_drhoa
2367
2368 exc_rhoa_norm_drhoa = scale_x*(e_lsda_x_arhoa*gx_anorm_drhoa + &
2369 e_lsda_x_a*(u_x_arhoanorm_drhoa*t18 + 0.2e1_dp*u_x_arhoa* &
2370 u_x_anorm_drhoa*c_x_2 + u_x_a*u_x_arhoanorm_drhoa*c_x_2)) + &
2371 scale_c*(e_lsda_c_abrhoa*gc_abnorm_drhoa + e_lsda_c_ab*( &
2372 u_c_abrhoanorm_drhoa*t170 + 0.2e1_dp*u_c_abrhoa* &
2373 u_c_abnorm_drhoa*c_cab_2 + u_c_ab*u_c_abrhoanorm_drhoa*c_cab_2 &
2374 ) + e_lsda_c_arhoa*gc_anorm_drhoa + e_lsda_c_a*( &
2375 u_c_arhoanorm_drhoa*t173 + 0.2e1_dp*u_c_arhoa*u_c_anorm_drhoa &
2376 *c_css_2 + u_c_a*u_c_arhoanorm_drhoa*c_css_2))
2377 e_r_ndr(ii) = e_r_ndr(ii) + 0.5_dp*0.5_dp*exc_rhoa_norm_drhoa
2378
2379 u_c_abrhobnorm_drhoa = -0.2e1_dp*t1436*t512 + 0.2e1_dp*t911 &
2380 *t1440*s_avg_2norm_drhoa
2381
2382 exc_rhob_norm_drhoa = scale_c*(e_lsda_c_abrhob*gc_abnorm_drhoa &
2383 + e_lsda_c_ab*(u_c_abrhobnorm_drhoa*t170 + 0.2e1_dp* &
2384 u_c_abrhob*u_c_abnorm_drhoa*c_cab_2 + u_c_ab* &
2385 u_c_abrhobnorm_drhoa*c_cab_2))
2386 e_r_ndr(ii) = e_r_ndr(ii) + 0.5_dp*0.5_dp*exc_rhob_norm_drhoa
2387
2388 t1571 = s_anorm_drhoa**2
2389 u_x_anorm_drhoanorm_drhoa = 0.2e1_dp*gamma_x*t1571*t16 - &
2390 0.10e2_dp*t568*t198*t1571 + 0.8e1_dp*t577*t579*t1571
2391 s_a_2norm_drhoanorm_drhoa = 0.2e1_dp*t1571
2392 s_a_21norm_drhoa = s_a_2norm_drhoa
2393 s_avg_2norm_drhoanorm_drhoa = s_a_2norm_drhoanorm_drhoa/0.2e1_dp
2394 s_avg_21norm_drhoa = s_a_21norm_drhoa/0.2e1_dp
2395 t1589 = t336*s_avg_2norm_drhoa
2396 t1590 = t339*s_avg_21norm_drhoa
2397 t1593 = t913*s_avg_2norm_drhoa
2398 u_c_abnorm_drhoanorm_drhoa = gamma_c_ab* &
2399 s_avg_2norm_drhoanorm_drhoa*t162 - 0.2e1_dp*t1589*t1590 + &
2400 0.2e1_dp*t911*t1593*s_avg_21norm_drhoa - t337*t339* &
2401 s_avg_2norm_drhoanorm_drhoa
2402 t1605 = t347*s_a_21norm_drhoa
2403 u_c_anorm_drhoanorm_drhoa = gamma_c_ss*s_a_2norm_drhoanorm_drhoa &
2404 *t165 - 0.2e1_dp*t344*s_a_2norm_drhoa*t1605 + 0.2e1_dp* &
2405 t930*t932*s_a_2norm_drhoa*s_a_21norm_drhoa - t345*t347* &
2406 s_a_2norm_drhoanorm_drhoa
2407
2408 exc_norm_drhoa_norm_drhoa = scale_x*e_lsda_x_a*( &
2409 u_x_anorm_drhoanorm_drhoa*t18 + 0.2e1_dp*u_x_anorm_drhoa**2* &
2410 c_x_2 + u_x_a*u_x_anorm_drhoanorm_drhoa*c_x_2) + scale_c*( &
2411 e_lsda_c_ab*(u_c_abnorm_drhoanorm_drhoa*t170 + 0.2e1_dp* &
2412 u_c_abnorm_drhoa*(gamma_c_ab*s_avg_21norm_drhoa*t162 - t337* &
2413 t1590)*c_cab_2 + u_c_ab*u_c_abnorm_drhoanorm_drhoa*c_cab_2) + &
2414 e_lsda_c_a*(u_c_anorm_drhoanorm_drhoa*t173 + 0.2e1_dp* &
2415 u_c_anorm_drhoa*(gamma_c_ss*s_a_21norm_drhoa*t165 - t345* &
2416 t1605)*c_css_2 + u_c_a*u_c_anorm_drhoanorm_drhoa*c_css_2))
2417 e_ndr_ndr(ii) = e_ndr_ndr(ii) + 0.5_dp*0.5_dp*exc_norm_drhoa_norm_drhoa
2418
2419 u_c_abrhoanorm_drhob = -0.2e1_dp*t906*t539 + 0.2e1_dp*t911* &
2420 t914*s_avg_2norm_drhob
2421
2422 exc_rhoa_norm_drhob = scale_c*(e_lsda_c_abrhoa*gc_abnorm_drhob &
2423 + e_lsda_c_ab*(u_c_abrhoanorm_drhob*t170 + 0.2e1_dp* &
2424 u_c_abrhoa*u_c_abnorm_drhob*c_cab_2 + u_c_ab* &
2425 u_c_abrhoanorm_drhob*c_cab_2))
2426 e_r_ndr(ii) = e_r_ndr(ii) + 0.5_dp*0.5_dp*exc_rhoa_norm_drhob
2427
2428 s_brhobnorm_drhob = -0.4e1_dp/0.3e1_dp*t367
2429 u_x_brhobnorm_drhob = 0.2e1_dp*gamma_x*s_bnorm_drhob*t371 - &
2430 0.10e2_dp*t1165*t377*s_bnorm_drhob + 0.2e1_dp*t370*t29* &
2431 s_brhobnorm_drhob + 0.8e1_dp*t1173*t1175*s_brhob* &
2432 s_bnorm_drhob - 0.2e1_dp*t374*t376*s_brhobnorm_drhob
2433 s_b_2rhobnorm_drhob = 0.2e1_dp*s_bnorm_drhob*s_brhob + &
2434 0.2e1_dp*s_b*s_brhobnorm_drhob
2435 s_avg_2rhobnorm_drhob = s_b_2rhobnorm_drhob/0.2e1_dp
2436 u_c_abrhobnorm_drhob = gamma_c_ab*s_avg_2rhobnorm_drhob*t162 - &
2437 0.2e1_dp*t1436*t539 + 0.2e1_dp*t911*t1440* &
2438 s_avg_2norm_drhob - t337*t339*s_avg_2rhobnorm_drhob
2439 u_c_brhobnorm_drhob = gamma_c_ss*s_b_2rhobnorm_drhob*t168 - &
2440 0.2e1_dp*t1451*t543 + 0.2e1_dp*t1455*t1458*s_b_2norm_drhob &
2441 - t484*t486*s_b_2rhobnorm_drhob
2442
2443 exc_rhob_norm_drhob = scale_x*(e_lsda_x_brhob*gx_bnorm_drhob + &
2444 e_lsda_x_b*(u_x_brhobnorm_drhob*t31 + 0.2e1_dp*u_x_brhob* &
2445 u_x_bnorm_drhob*c_x_2 + u_x_b*u_x_brhobnorm_drhob*c_x_2)) + &
2446 scale_c*(e_lsda_c_abrhob*gc_abnorm_drhob + e_lsda_c_ab*( &
2447 u_c_abrhobnorm_drhob*t170 + 0.2e1_dp*u_c_abrhob* &
2448 u_c_abnorm_drhob*c_cab_2 + u_c_ab*u_c_abrhobnorm_drhob*c_cab_2 &
2449 ) + e_lsda_c_brhob*gc_bnorm_drhob + e_lsda_c_b*( &
2450 u_c_brhobnorm_drhob*t176 + 0.2e1_dp*u_c_brhob*u_c_bnorm_drhob &
2451 *c_css_2 + u_c_b*u_c_brhobnorm_drhob*c_css_2))
2452 e_r_ndr(ii) = e_r_ndr(ii) + 0.5_dp*0.5_dp*exc_rhob_norm_drhob
2453
2454 u_c_abnorm_drhoanorm_drhob = -0.2e1_dp*t1589*t539 + 0.2e1_dp* &
2455 t911*t1593*s_avg_2norm_drhob
2456
2457 exc_norm_drhoa_norm_drhob = scale_c*e_lsda_c_ab*( &
2458 u_c_abnorm_drhoanorm_drhob*t170 + 0.2e1_dp*u_c_abnorm_drhoa* &
2459 u_c_abnorm_drhob*c_cab_2 + u_c_ab*u_c_abnorm_drhoanorm_drhob* &
2460 c_cab_2)
2461 e_ndr_ndr(ii) = e_ndr_ndr(ii) + 0.5_dp*exc_norm_drhoa_norm_drhob
2462
2463 t1719 = s_bnorm_drhob**2
2464 u_x_bnorm_drhobnorm_drhob = 0.2e1_dp*gamma_x*t1719*t29 - &
2465 0.10e2_dp*t1165*t376*t1719 + 0.8e1_dp*t1173*t1175*t1719
2466 s_b_2norm_drhobnorm_drhob = 0.2e1_dp*t1719
2467 s_b_21norm_drhob = s_b_2norm_drhob
2468 s_avg_2norm_drhobnorm_drhob = s_b_2norm_drhobnorm_drhob/0.2e1_dp
2469 s_avg_21norm_drhob = s_b_21norm_drhob/0.2e1_dp
2470 t1738 = t339*s_avg_21norm_drhob
2471 u_c_abnorm_drhobnorm_drhob = gamma_c_ab* &
2472 s_avg_2norm_drhobnorm_drhob*t162 - 0.2e1_dp*t336* &
2473 s_avg_2norm_drhob*t1738 + 0.2e1_dp*t911*t913* &
2474 s_avg_2norm_drhob*s_avg_21norm_drhob - t337*t339* &
2475 s_avg_2norm_drhobnorm_drhob
2476 t1753 = t486*s_b_21norm_drhob
2477 u_c_bnorm_drhobnorm_drhob = gamma_c_ss*s_b_2norm_drhobnorm_drhob &
2478 *t168 - 0.2e1_dp*t344*s_b_2norm_drhob*t1753 + 0.2e1_dp* &
2479 t1455*t1457*s_b_2norm_drhob*s_b_21norm_drhob - t484*t486* &
2480 s_b_2norm_drhobnorm_drhob
2481
2482 exc_norm_drhob_norm_drhob = scale_x*e_lsda_x_b*( &
2483 u_x_bnorm_drhobnorm_drhob*t31 + 0.2e1_dp*u_x_bnorm_drhob**2* &
2484 c_x_2 + u_x_b*u_x_bnorm_drhobnorm_drhob*c_x_2) + scale_c*( &
2485 e_lsda_c_ab*(u_c_abnorm_drhobnorm_drhob*t170 + 0.2e1_dp* &
2486 u_c_abnorm_drhob*(gamma_c_ab*s_avg_21norm_drhob*t162 - t337* &
2487 t1738)*c_cab_2 + u_c_ab*u_c_abnorm_drhobnorm_drhob*c_cab_2) + &
2488 e_lsda_c_b*(u_c_bnorm_drhobnorm_drhob*t176 + 0.2e1_dp* &
2489 u_c_bnorm_drhob*(gamma_c_ss*s_b_21norm_drhob*t168 - t484* &
2490 t1753)*c_css_2 + u_c_b*u_c_bnorm_drhobnorm_drhob*c_css_2))
2491 e_ndr_ndr(ii) = e_ndr_ndr(ii) + 0.5_dp*0.5_dp*exc_norm_drhob_norm_drhob
2492 END IF ! <1 || >1
2493 END IF ! /=0
2494 END IF ! rho>epsilon_rho
2495 END DO
2496!$OMP END DO
2497 END SELECT
2498 END SUBROUTINE b97_lda_calc
2499
2500END MODULE xc_b97
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::grimme2006
integer, save, public grimme2006
Definition bibliography.F:36
bibliography::becke1997
integer, save, public becke1997
Definition bibliography.F:36
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
xc_b97
calculates the b97 correlation functional
Definition xc_b97.F:24
xc_b97::b97_lda_info
subroutine, public b97_lda_info(b97_params, reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_b97.F:187
xc_b97::b97_lda_eval
subroutine, public b97_lda_eval(rho_set, deriv_set, grad_deriv, b97_params)
evaluates the b97 correlation functional for lda
Definition xc_b97.F:250
xc_b97::b97_lsd_info
subroutine, public b97_lsd_info(b97_params, reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_b97.F:219
xc_b97::b97_lsd_eval
subroutine, public b97_lsd_eval(rho_set, deriv_set, grad_deriv, b97_params)
evaluates the b 97 xc functional for lsd
Definition xc_b97.F:360
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition xc_derivative_desc.F:26
xc_derivative_desc::deriv_norm_drho
integer, parameter, public deriv_norm_drho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
Definition xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:65
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::xc_b97_mardirossian
integer, parameter, public xc_b97_mardirossian
Definition xc_input_constants.F:96
xc_input_constants::xc_b97_orig
integer, parameter, public xc_b97_orig
Definition xc_input_constants.F:96
xc_input_constants::xc_b97_grimme
integer, parameter, public xc_b97_grimme
Definition xc_input_constants.F:96
xc_input_constants::xc_b97_3c
integer, parameter, public xc_b97_3c
Definition xc_input_constants.F:96
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition xc_rho_set_types.F:281
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78