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