Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief calculates the 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 623 : 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 623 : IF (lda) THEN
88 623 : pol = "(unpolarized)"
89 : ELSE
90 0 : pol = "(polarized)"
91 : END IF
92 623 : SELECT CASE (param)
93 : CASE (xc_b97_orig)
94 0 : CALL cite_reference(Becke1997)
95 0 : IF (sx == 1._dp .AND. sc == 1._dp) THEN
96 0 : IF (PRESENT(reference)) THEN
97 : reference = "A.D.Becke, "// &
98 : "J. Chem. Phys, vol. 107, pp. 8554, (1997), needs exact x, "// &
99 0 : pol
100 : END IF
101 0 : IF (PRESENT(shortform)) THEN
102 0 : shortform = "B97, Becke 1997 xc functional, needs exact x "//pol
103 : END IF
104 : ELSE
105 0 : IF (PRESENT(reference)) THEN
106 : WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a,a)") &
107 0 : "A.D.Becke, ", &
108 0 : "J. Chem. Phys, vol. 107, pp. 8554, (1997)", &
109 0 : sx, sc, ", needs exact x ", pol
110 : END IF
111 0 : IF (PRESENT(shortform)) THEN
112 : WRITE (shortform, "(a,a,'sx=',f5.3,'sc=',f5.3)") &
113 0 : "B97, Becke 1997 xc functional (unpolarized)", pol, sx, sc
114 : END IF
115 : END IF
116 : CASE (xc_b97_grimme)
117 579 : CALL cite_reference(Becke1997)
118 579 : CALL cite_reference(Grimme2006)
119 579 : IF (sx == 1._dp .AND. sc == 1._dp) THEN
120 0 : IF (PRESENT(reference)) THEN
121 : reference = "B97-D, Grimme xc functional,"// &
122 : " J Comput Chem 27: 1787-1799 (2006),"// &
123 0 : " needs C6 dispersion "//pol
124 : END IF
125 0 : IF (PRESENT(shortform)) THEN
126 0 : shortform = "B97-D, Grimme xc functional "//pol
127 : END IF
128 : ELSE
129 579 : IF (PRESENT(reference)) THEN
130 : WRITE (reference, "(a,a,3x,' sx=',f6.3,' sc=',f6.3,1x,a)") &
131 4 : "B97-D, Grimme xc functional,", &
132 4 : " J Comput Chem 27: 1787-1799 (2006) ", &
133 8 : sx, sc, pol
134 : END IF
135 579 : IF (PRESENT(shortform)) THEN
136 : WRITE (shortform, "(a,a,3x,' sx=',f6.3,' sc=',f6.3,' (LDA)')") &
137 4 : "B97-D, Grimme xc functional ", pol, sx, sc
138 : END IF
139 : END IF
140 : CASE (xc_b97_mardirossian)
141 0 : CALL cite_reference(Becke1997)
142 : ! CALL cite_reference(Mardirossian2014)
143 0 : IF (sx == 1._dp .AND. sc == 1._dp) THEN
144 0 : IF (PRESENT(reference)) THEN
145 : reference = "wB97X-V, xc functional,"// &
146 : " Mardirossian and Head-Gordon, PCCP DOI: 10.1039/c3cp54374a,"// &
147 0 : " needs HFX exchange and VV10 dispersion (NOT TESTED)"//pol
148 : END IF
149 0 : IF (PRESENT(shortform)) THEN
150 0 : shortform = "wB97X-V, HFX+B97+VV10 functional (NOT TESTED)"//pol
151 : END IF
152 : ELSE
153 0 : CPABORT("Unsupported scaling factors")
154 : END IF
155 : CASE (xc_b97_3c)
156 44 : CALL cite_reference(Becke1997)
157 : ! CALL cite_reference(Mardirossian2014)
158 44 : IF (sx == 1._dp .AND. sc == 1._dp) THEN
159 44 : 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 0 : " with CALCULATE_C9_TERM and SRB correction"//pol
164 : END IF
165 44 : IF (PRESENT(shortform)) THEN
166 0 : shortform = "B97-3c composite method"//pol
167 : END IF
168 : ELSE
169 0 : CPABORT("Unsupported scaling factors")
170 : END IF
171 : CASE default
172 623 : CPABORT("Unsupported parametrization")
173 : END SELECT
174 623 : 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 1869 : 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 623 : CALL section_vals_val_get(b97_params, "parametrization", i_val=param)
196 623 : CALL section_vals_val_get(b97_params, "scale_x", r_val=sx)
197 623 : CALL section_vals_val_get(b97_params, "scale_c", r_val=sc)
198 :
199 3099 : CALL b97_ref(param, .TRUE., sx, sc, reference, shortform)
200 623 : IF (PRESENT(needs)) THEN
201 619 : needs%rho = .TRUE.
202 619 : needs%norm_drho = .TRUE.
203 : END IF
204 623 : IF (PRESENT(max_deriv)) max_deriv = 2
205 :
206 623 : 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 0 : 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 0 : CALL section_vals_val_get(b97_params, "parametrization", i_val=param)
228 0 : CALL section_vals_val_get(b97_params, "scale_x", r_val=sx)
229 0 : CALL section_vals_val_get(b97_params, "scale_c", r_val=sc)
230 :
231 0 : CALL b97_ref(param, .FALSE., sx, sc, reference, shortform)
232 0 : IF (PRESENT(needs)) THEN
233 0 : needs%rho_spin = .TRUE.
234 0 : needs%norm_drho_spin = .TRUE.
235 : END IF
236 0 : IF (PRESENT(max_deriv)) max_deriv = 2
237 :
238 0 : 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 3035 : 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 607 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
264 607 : e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
265 607 : e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
266 : TYPE(xc_derivative_type), POINTER :: deriv
267 :
268 607 : 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 607 : drho_cutoff=epsilon_norm_drho)
273 607 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
274 :
275 607 : dummy => rho
276 :
277 607 : e_0 => dummy
278 607 : e_rho => dummy
279 607 : e_ndrho => dummy
280 607 : e_rho_rho => dummy
281 607 : e_ndrho_rho => dummy
282 607 : e_ndrho_ndrho => dummy
283 607 : e_rho_rho_rho => dummy
284 607 : e_ndrho_rho_rho => dummy
285 607 : e_ndrho_ndrho_rho => dummy
286 607 : e_ndrho_ndrho_ndrho => dummy
287 :
288 607 : IF (grad_deriv >= 0) THEN
289 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
290 607 : allocate_deriv=.TRUE.)
291 607 : CALL xc_derivative_get(deriv, deriv_data=e_0)
292 : END IF
293 607 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
294 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
295 607 : allocate_deriv=.TRUE.)
296 607 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
297 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
298 607 : allocate_deriv=.TRUE.)
299 607 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
300 : END IF
301 607 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
302 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
303 0 : allocate_deriv=.TRUE.)
304 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
305 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
306 0 : allocate_deriv=.TRUE.)
307 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
308 : deriv => xc_dset_get_derivative(deriv_set, &
309 0 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
310 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
311 : END IF
312 607 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
313 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
314 0 : allocate_deriv=.TRUE.)
315 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
316 : deriv => xc_dset_get_derivative(deriv_set, &
317 0 : [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.TRUE.)
318 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
319 : deriv => xc_dset_get_derivative(deriv_set, &
320 0 : [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.TRUE.)
321 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
322 : deriv => xc_dset_get_derivative(deriv_set, &
323 0 : [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
324 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
325 : END IF
326 607 : IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
327 0 : CPABORT("derivatives bigger than 3 not implemented")
328 : END IF
329 :
330 607 : CALL section_vals_val_get(b97_params, "parametrization", i_val=param)
331 607 : CALL section_vals_val_get(b97_params, "scale_c", r_val=scale_c)
332 607 : 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 607 : !$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 607 : NULLIFY (dummy)
352 607 : CALL timestop(handle)
353 607 : 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 0 : 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 0 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndra, e_ndra_ndra, &
375 0 : e_ndra_ndrb, e_ndra_ra, e_ndra_rb, e_ndrb, e_ndrb_ndrb, e_ndrb_ra, e_ndrb_rb, e_ra, &
376 0 : 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 0 : CALL timeset(routineN, handle)
380 0 : 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 0 : local_bounds=bo)
387 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
388 :
389 0 : dummy => rhoa
390 0 : e_0 => dummy
391 0 : e_ra => dummy
392 0 : e_rb => dummy
393 0 : e_ndra => dummy
394 0 : e_ndrb => dummy
395 :
396 0 : e_ndra_ra => dummy
397 0 : e_ndra_rb => dummy
398 0 : e_ndrb_rb => dummy
399 0 : e_ndrb_ra => dummy
400 :
401 0 : e_ndra_ndra => dummy
402 0 : e_ndrb_ndrb => dummy
403 0 : e_ndra_ndrb => dummy
404 :
405 0 : e_ra_ra => dummy
406 0 : e_ra_rb => dummy
407 0 : e_rb_rb => dummy
408 :
409 0 : IF (grad_deriv >= 0) THEN
410 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
411 0 : allocate_deriv=.TRUE.)
412 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
413 : END IF
414 0 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
415 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
416 0 : allocate_deriv=.TRUE.)
417 0 : CALL xc_derivative_get(deriv, deriv_data=e_ra)
418 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
419 0 : allocate_deriv=.TRUE.)
420 0 : CALL xc_derivative_get(deriv, deriv_data=e_rb)
421 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
422 0 : allocate_deriv=.TRUE.)
423 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndra)
424 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
425 0 : allocate_deriv=.TRUE.)
426 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
427 : END IF
428 0 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
429 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
430 0 : allocate_deriv=.TRUE.)
431 0 : CALL xc_derivative_get(deriv, deriv_data=e_ra_ra)
432 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
433 0 : allocate_deriv=.TRUE.)
434 0 : CALL xc_derivative_get(deriv, deriv_data=e_ra_rb)
435 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
436 0 : allocate_deriv=.TRUE.)
437 0 : CALL xc_derivative_get(deriv, deriv_data=e_rb_rb)
438 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
439 0 : allocate_deriv=.TRUE.)
440 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_ra)
441 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
442 0 : allocate_deriv=.TRUE.)
443 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_rb)
444 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
445 0 : allocate_deriv=.TRUE.)
446 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb_rb)
447 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa], &
448 0 : allocate_deriv=.TRUE.)
449 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb_ra)
450 : deriv => xc_dset_get_derivative(deriv_set, &
451 0 : [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
452 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndra)
453 : deriv => xc_dset_get_derivative(deriv_set, &
454 0 : [deriv_norm_drhoa, deriv_norm_drhob], allocate_deriv=.TRUE.)
455 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndrb)
456 : deriv => xc_dset_get_derivative(deriv_set, &
457 0 : [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
458 0 : 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 0 : CALL section_vals_val_get(b97_params, "parametrization", i_val=param)
465 0 : CALL section_vals_val_get(b97_params, "scale_x", r_val=scale_x)
466 0 : 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 0 : !$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 0 : NULLIFY (dummy)
493 0 : CALL timestop(handle)
494 0 : END SUBROUTINE b97_lsd_eval
495 :
496 : ! **************************************************************************************************
497 : !> \brief ...
498 : !> \param param ...
499 : !> \return ...
500 : ! **************************************************************************************************
501 607 : FUNCTION b97_coeffs(param) RESULT(res)
502 : INTEGER, INTENT(in) :: param
503 : REAL(dp), DIMENSION(10) :: res
504 :
505 607 : SELECT CASE (param)
506 : CASE (xc_b97_orig)
507 0 : res = params_b97_orig
508 : CASE (xc_b97_grimme)
509 6281 : res = params_b97_grimme
510 : CASE (xc_b97_mardirossian)
511 0 : res = params_b97_mardirossian
512 : CASE (xc_b97_3c)
513 396 : res = params_b97_3c
514 : CASE default
515 0 : CPABORT("Unsupported parametrization")
516 607 : res = 0.0_dp
517 : END SELECT
518 607 : 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 0 : 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 0 : p_1 = 0.10e1_dp
633 0 : A_1 = 0.31091e-1_dp
634 0 : alpha_1_1 = 0.21370e0_dp
635 0 : beta_1_1 = 0.75957e1_dp
636 0 : beta_2_1 = 0.35876e1_dp
637 0 : beta_3_1 = 0.16382e1_dp
638 0 : beta_4_1 = 0.49294e0_dp
639 0 : p_2 = 0.10e1_dp
640 0 : A_2 = 0.15545e-1_dp
641 0 : alpha_1_2 = 0.20548e0_dp
642 0 : beta_1_2 = 0.141189e2_dp
643 0 : beta_2_2 = 0.61977e1_dp
644 0 : beta_3_2 = 0.33662e1_dp
645 0 : beta_4_2 = 0.62517e0_dp
646 0 : p_3 = 0.10e1_dp
647 0 : A_3 = 0.16887e-1_dp
648 0 : alpha_1_3 = 0.11125e0_dp
649 0 : beta_1_3 = 0.10357e2_dp
650 0 : beta_2_3 = 0.36231e1_dp
651 0 : beta_3_3 = 0.88026e0_dp
652 0 : beta_4_3 = 0.49671e0_dp
653 :
654 0 : coeffs = b97_coeffs(param)
655 0 : c_x_0 = coeffs(1)
656 0 : c_x_1 = coeffs(2)
657 0 : c_x_2 = coeffs(3)
658 0 : c_cab_0 = coeffs(4)
659 0 : c_cab_1 = coeffs(5)
660 0 : c_cab_2 = coeffs(6)
661 0 : c_css_0 = coeffs(7)
662 0 : c_css_1 = coeffs(8)
663 0 : c_css_2 = coeffs(9)
664 :
665 0 : scale_c = scale_c_in
666 0 : scale_x = scale_x_in
667 0 : IF (scale_x == -1.0_dp) scale_x = coeffs(10)
668 :
669 0 : gamma_x = 0.004_dp
670 0 : gamma_c_ss = 0.2_dp
671 0 : gamma_c_ab = 0.006_dp
672 :
673 0 : SELECT CASE (grad_deriv)
674 : CASE default
675 0 : t1 = 3**(0.1e1_dp/0.3e1_dp)
676 0 : t2 = 4**(0.1e1_dp/0.3e1_dp)
677 0 : t3 = t2**2
678 0 : t4 = t1*t3
679 0 : t6 = (0.1e1_dp/pi)**(0.1e1_dp/0.3e1_dp)
680 0 : !$OMP DO
681 : DO ii = 1, npoints
682 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
683 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
684 0 : rho = my_rhoa + my_rhob
685 0 : IF (rho > epsilon_rho) THEN
686 0 : my_rhoa = MAX(my_rhoa, 0.5_dp*epsilon_rho)
687 0 : my_rhob = MAX(my_rhob, 0.5_dp*epsilon_rho)
688 0 : rho = my_rhoa + my_rhob
689 0 : my_norm_drhoa = norm_drhoa(ii)
690 0 : my_norm_drhob = norm_drhob(ii)
691 0 : t7 = my_rhoa**(0.1e1_dp/0.3e1_dp)
692 0 : t8 = t7*my_rhoa
693 0 : e_lsda_x_a = -0.3e1_dp/0.8e1_dp*t4*t6*t8
694 0 : t12 = 0.1e1_dp/t8
695 0 : s_a = my_norm_drhoa*t12
696 0 : t13 = s_a**2
697 0 : t14 = gamma_x*t13
698 0 : t15 = 0.1e1_dp + t14
699 0 : t16 = 0.1e1_dp/t15
700 0 : u_x_a = t14*t16
701 0 : t18 = c_x_1 + u_x_a*c_x_2
702 0 : gx_a = c_x_0 + u_x_a*t18
703 0 : t20 = my_rhob**(0.1e1_dp/0.3e1_dp)
704 0 : t21 = t20*my_rhob
705 0 : e_lsda_x_b = -0.3e1_dp/0.8e1_dp*t4*t6*t21
706 0 : t25 = 0.1e1_dp/t21
707 0 : s_b = my_norm_drhob*t25
708 0 : t26 = s_b**2
709 0 : t27 = gamma_x*t26
710 0 : t28 = 0.1e1_dp + t27
711 0 : t29 = 0.1e1_dp/t28
712 0 : u_x_b = t27*t29
713 0 : t31 = c_x_1 + u_x_b*c_x_2
714 0 : gx_b = c_x_0 + u_x_b*t31
715 0 : t33 = my_rhoa - my_rhob
716 0 : t34 = 0.1e1_dp/rho
717 0 : chi = t33*t34
718 0 : t35 = 0.1e1_dp/pi
719 0 : t36 = t35*t34
720 0 : t37 = t36**(0.1e1_dp/0.3e1_dp)
721 0 : rs = t4*t37/0.4e1_dp
722 0 : t40 = 0.1e1_dp + alpha_1_1*rs
723 0 : t42 = 0.1e1_dp/A_1
724 0 : t43 = SQRT(rs)
725 0 : t46 = t43*rs
726 0 : t48 = p_1 + 0.1e1_dp
727 0 : t49 = rs**t48
728 0 : t50 = beta_4_1*t49
729 0 : t51 = beta_1_1*t43 + beta_2_1*rs + beta_3_1*t46 + t50
730 0 : t55 = 0.1e1_dp + t42/t51/0.2e1_dp
731 0 : t56 = LOG(t55)
732 0 : e_c_u_0 = -0.2e1_dp*A_1*t40*t56
733 0 : t60 = 0.1e1_dp + alpha_1_2*rs
734 0 : t62 = 0.1e1_dp/A_2
735 0 : t66 = p_2 + 0.1e1_dp
736 0 : t67 = rs**t66
737 0 : t68 = beta_4_2*t67
738 0 : t69 = beta_1_2*t43 + beta_2_2*rs + beta_3_2*t46 + t68
739 0 : t73 = 0.1e1_dp + t62/t69/0.2e1_dp
740 0 : t74 = LOG(t73)
741 0 : t78 = 0.1e1_dp + alpha_1_3*rs
742 0 : t80 = 0.1e1_dp/A_3
743 0 : t84 = p_3 + 0.1e1_dp
744 0 : t85 = rs**t84
745 0 : t86 = beta_4_3*t85
746 0 : t87 = beta_1_3*t43 + beta_2_3*rs + beta_3_3*t46 + t86
747 0 : t91 = 0.1e1_dp + t80/t87/0.2e1_dp
748 0 : t92 = LOG(t91)
749 0 : alpha_c = 0.2e1_dp*A_3*t78*t92
750 0 : t94 = 2**(0.1e1_dp/0.3e1_dp)
751 0 : t97 = 1/(2*t94 - 2)
752 0 : t98 = 0.1e1_dp + chi
753 0 : t99 = t98**(0.1e1_dp/0.3e1_dp)
754 0 : t101 = 0.1e1_dp - chi
755 0 : t102 = t101**(0.1e1_dp/0.3e1_dp)
756 0 : f = (t99*t98 + t102*t101 - 0.2e1_dp)*t97
757 0 : t105 = alpha_c*f
758 0 : t106 = 0.9e1_dp/0.8e1_dp/t97
759 0 : t107 = chi**2
760 0 : t108 = t107**2
761 0 : t110 = t106*(0.1e1_dp - t108)
762 0 : t112 = -0.2e1_dp*A_2*t60*t74 - e_c_u_0
763 0 : t113 = t112*f
764 0 : epsilon_c_unif = e_c_u_0 + t105*t110 + t113*t108
765 0 : t116 = t35/my_rhoa
766 0 : t117 = t116**(0.1e1_dp/0.3e1_dp)
767 0 : rs_a = t4*t117/0.4e1_dp
768 0 : t120 = 0.1e1_dp + alpha_1_2*rs_a
769 0 : t122 = SQRT(rs_a)
770 0 : t125 = t122*rs_a
771 0 : t127 = rs_a**t66
772 0 : t128 = beta_4_2*t127
773 0 : t129 = beta_1_2*t122 + beta_2_2*rs_a + beta_3_2*t125 + t128
774 0 : t133 = 0.1e1_dp + t62/t129/0.2e1_dp
775 0 : t134 = LOG(t133)
776 0 : epsilon_c_unif_a = -0.2e1_dp*A_2*t120*t134
777 0 : t138 = t35/my_rhob
778 0 : t139 = t138**(0.1e1_dp/0.3e1_dp)
779 0 : rs_b = t4*t139/0.4e1_dp
780 0 : t142 = 0.1e1_dp + alpha_1_2*rs_b
781 0 : t144 = SQRT(rs_b)
782 0 : t147 = t144*rs_b
783 0 : t149 = rs_b**t66
784 0 : t150 = beta_4_2*t149
785 0 : t151 = beta_1_2*t144 + beta_2_2*rs_b + beta_3_2*t147 + t150
786 0 : t155 = 0.1e1_dp + t62/t151/0.2e1_dp
787 0 : t156 = LOG(t155)
788 0 : epsilon_c_unif_b = -0.2e1_dp*A_2*t142*t156
789 0 : s_a_2 = t13
790 0 : s_b_2 = t26
791 0 : s_avg_2 = s_a_2/0.2e1_dp + s_b_2/0.2e1_dp
792 0 : e_lsda_c_a = epsilon_c_unif_a*my_rhoa
793 0 : e_lsda_c_b = epsilon_c_unif_b*my_rhob
794 0 : t160 = gamma_c_ab*s_avg_2
795 0 : t161 = 0.1e1_dp + t160
796 0 : t162 = 0.1e1_dp/t161
797 0 : u_c_ab = t160*t162
798 0 : t163 = gamma_c_ss*s_a_2
799 0 : t164 = 0.1e1_dp + t163
800 0 : t165 = 0.1e1_dp/t164
801 0 : u_c_a = t163*t165
802 0 : t166 = gamma_c_ss*s_b_2
803 0 : t167 = 0.1e1_dp + t166
804 0 : t168 = 0.1e1_dp/t167
805 0 : u_c_b = t166*t168
806 0 : e_lsda_c_ab = epsilon_c_unif*rho - e_lsda_c_a - e_lsda_c_b
807 0 : t170 = c_cab_1 + u_c_ab*c_cab_2
808 0 : gc_ab = c_cab_0 + u_c_ab*t170
809 0 : t173 = c_css_1 + u_c_a*c_css_2
810 0 : gc_a = c_css_0 + u_c_a*t173
811 0 : t176 = c_css_1 + u_c_b*c_css_2
812 0 : gc_b = c_css_0 + u_c_b*t176
813 :
814 0 : IF (grad_deriv >= 0) THEN
815 : exc = scale_x*(e_lsda_x_a*gx_a + e_lsda_x_b*gx_b) + scale_c &
816 0 : *(e_lsda_c_ab*gc_ab + e_lsda_c_a*gc_a + e_lsda_c_b*gc_b)
817 0 : e_0(ii) = e_0(ii) + exc
818 : END IF
819 :
820 0 : IF (grad_deriv /= 0) THEN
821 :
822 0 : e_lsda_x_arhoa = -t4*t6*t7/0.2e1_dp
823 0 : t186 = my_rhoa**2
824 0 : t188 = 0.1e1_dp/t7/t186
825 0 : s_arhoa = -0.4e1_dp/0.3e1_dp*my_norm_drhoa*t188
826 0 : t191 = gamma_x*s_a
827 0 : t192 = t16*s_arhoa
828 0 : t194 = gamma_x**2
829 0 : t196 = t194*s_a_2*s_a
830 0 : t197 = t15**2
831 0 : t198 = 0.1e1_dp/t197
832 0 : t199 = t198*s_arhoa
833 0 : u_x_arhoa = 0.2e1_dp*t191*t192 - 0.2e1_dp*t196*t199
834 0 : gx_arhoa = u_x_arhoa*t18 + u_x_a*u_x_arhoa*c_x_2
835 0 : t207 = rho**2
836 0 : t208 = 0.1e1_dp/t207
837 0 : t209 = t33*t208
838 0 : chirhoa = t34 - t209
839 0 : t210 = t37**2
840 0 : t212 = 0.1e1_dp/t210*t35
841 0 : rsrhoa = -t4*t212*t208/0.12e2_dp
842 0 : t216 = A_1*alpha_1_1
843 0 : t220 = t51**2
844 0 : t221 = 0.1e1_dp/t220
845 0 : t222 = t40*t221
846 0 : t223 = 0.1e1_dp/t43
847 0 : t224 = beta_1_1*t223
848 0 : t228 = beta_3_1*t43
849 0 : t232 = 0.1e1_dp/rs
850 : t235 = t224*rsrhoa/0.2e1_dp + beta_2_1*rsrhoa + 0.3e1_dp/ &
851 0 : 0.2e1_dp*t228*rsrhoa + t50*t48*rsrhoa*t232
852 0 : t236 = 0.1e1_dp/t55
853 0 : t237 = t235*t236
854 0 : e_c_u_0rhoa = -0.2e1_dp*t216*rsrhoa*t56 + t222*t237
855 0 : t239 = A_2*alpha_1_2
856 0 : t243 = t69**2
857 0 : t244 = 0.1e1_dp/t243
858 0 : t245 = t60*t244
859 0 : t246 = beta_1_2*t223
860 0 : t250 = beta_3_2*t43
861 : t256 = t246*rsrhoa/0.2e1_dp + beta_2_2*rsrhoa + 0.3e1_dp/ &
862 0 : 0.2e1_dp*t250*rsrhoa + t68*t66*rsrhoa*t232
863 0 : t257 = 0.1e1_dp/t73
864 0 : t258 = t256*t257
865 0 : e_c_u_1rhoa = -0.2e1_dp*t239*rsrhoa*t74 + t245*t258
866 0 : t260 = A_3*alpha_1_3
867 0 : t264 = t87**2
868 0 : t265 = 0.1e1_dp/t264
869 0 : t266 = t78*t265
870 0 : t267 = beta_1_3*t223
871 0 : t271 = beta_3_3*t43
872 : t277 = t267*rsrhoa/0.2e1_dp + beta_2_3*rsrhoa + 0.3e1_dp/ &
873 0 : 0.2e1_dp*t271*rsrhoa + t86*t84*rsrhoa*t232
874 0 : t278 = 0.1e1_dp/t91
875 0 : t279 = t277*t278
876 0 : 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 0 : *t102*chirhoa)*t97
879 0 : t285 = alpha_crhoa*f
880 0 : t287 = alpha_c*frhoa
881 0 : t289 = t107*chi
882 0 : t290 = t106*t289
883 0 : t291 = t290*chirhoa
884 0 : t293 = 0.4e1_dp*t105*t291
885 0 : t294 = e_c_u_1rhoa - e_c_u_0rhoa
886 0 : t295 = t294*f
887 0 : t297 = t112*frhoa
888 0 : t299 = t289*chirhoa
889 0 : t301 = 0.4e1_dp*t113*t299
890 : epsilon_c_unifrhoa = e_c_u_0rhoa + t285*t110 + t287*t110 - &
891 0 : t293 + t295*t108 + t297*t108 + t301
892 0 : t302 = t117**2
893 0 : t304 = 0.1e1_dp/t302*t35
894 0 : rs_arhoa = -t4*t304/t186/0.12e2_dp
895 0 : t312 = t129**2
896 0 : t313 = 0.1e1_dp/t312
897 0 : t314 = t120*t313
898 0 : t315 = 0.1e1_dp/t122
899 0 : t316 = beta_1_2*t315
900 0 : t320 = beta_3_2*t122
901 0 : t324 = 0.1e1_dp/rs_a
902 : t327 = t316*rs_arhoa/0.2e1_dp + beta_2_2*rs_arhoa + 0.3e1_dp &
903 0 : /0.2e1_dp*t320*rs_arhoa + t128*t66*rs_arhoa*t324
904 0 : t328 = 0.1e1_dp/t133
905 : epsilon_c_unif_arhoa = -0.2e1_dp*t239*rs_arhoa*t134 + t314* &
906 0 : t327*t328
907 0 : s_a_2rhoa = 0.2e1_dp*s_a*s_arhoa
908 0 : s_avg_2rhoa = s_a_2rhoa/0.2e1_dp
909 0 : e_lsda_c_arhoa = epsilon_c_unif_arhoa*my_rhoa + epsilon_c_unif_a
910 0 : t336 = gamma_c_ab**2
911 0 : t337 = t336*s_avg_2
912 0 : t338 = t161**2
913 0 : t339 = 0.1e1_dp/t338
914 0 : u_c_abrhoa = gamma_c_ab*s_avg_2rhoa*t162 - t337*t339*s_avg_2rhoa
915 0 : t344 = gamma_c_ss**2
916 0 : t345 = t344*s_a_2
917 0 : t346 = t164**2
918 0 : t347 = 0.1e1_dp/t346
919 0 : 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 0 : e_lsda_c_arhoa
922 0 : gc_abrhoa = u_c_abrhoa*t170 + u_c_ab*u_c_abrhoa*c_cab_2
923 0 : gc_arhoa = u_c_arhoa*t173 + u_c_a*u_c_arhoa*c_css_2
924 :
925 0 : 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 0 : gc_abrhoa + e_lsda_c_arhoa*gc_a + e_lsda_c_a*gc_arhoa)
929 0 : e_ra(ii) = e_ra(ii) + exc_rhoa
930 : END IF
931 :
932 0 : e_lsda_x_brhob = -t4*t6*t20/0.2e1_dp
933 0 : t365 = my_rhob**2
934 0 : t367 = 0.1e1_dp/t20/t365
935 0 : s_brhob = -0.4e1_dp/0.3e1_dp*my_norm_drhob*t367
936 0 : t370 = gamma_x*s_b
937 0 : t371 = t29*s_brhob
938 0 : t374 = t194*s_b_2*s_b
939 0 : t375 = t28**2
940 0 : t376 = 0.1e1_dp/t375
941 0 : t377 = t376*s_brhob
942 0 : u_x_brhob = 0.2e1_dp*t370*t371 - 0.2e1_dp*t374*t377
943 0 : gx_brhob = u_x_brhob*t31 + u_x_b*u_x_brhob*c_x_2
944 0 : chirhob = -t34 - t209
945 0 : rsrhob = rsrhoa
946 : t396 = t224*rsrhob/0.2e1_dp + beta_2_1*rsrhob + 0.3e1_dp/ &
947 0 : 0.2e1_dp*t228*rsrhob + t50*t48*rsrhob*t232
948 0 : 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 : 0.2e1_dp*t250*rsrhob + t68*t66*rsrhob*t232
951 0 : 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 : 0.2e1_dp*t271*rsrhob + t86*t84*rsrhob*t232
954 0 : 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 0 : *t102*chirhob)*t97
957 0 : t431 = alpha_crhob*f
958 0 : t433 = alpha_c*frhob
959 0 : t435 = t290*chirhob
960 0 : t437 = 0.4e1_dp*t105*t435
961 0 : t438 = e_c_u_1rhob - e_c_u_0rhob
962 0 : t439 = t438*f
963 0 : t441 = t112*frhob
964 0 : t443 = t289*chirhob
965 0 : t445 = 0.4e1_dp*t113*t443
966 : epsilon_c_unifrhob = e_c_u_0rhob + t431*t110 + t433*t110 - &
967 0 : t437 + t439*t108 + t441*t108 + t445
968 0 : t446 = t139**2
969 0 : t448 = 0.1e1_dp/t446*t35
970 0 : rs_brhob = -t4*t448/t365/0.12e2_dp
971 0 : t456 = t151**2
972 0 : t457 = 0.1e1_dp/t456
973 0 : t458 = t142*t457
974 0 : t459 = 0.1e1_dp/t144
975 0 : t460 = beta_1_2*t459
976 0 : t464 = beta_3_2*t144
977 0 : t468 = 0.1e1_dp/rs_b
978 : t471 = t460*rs_brhob/0.2e1_dp + beta_2_2*rs_brhob + 0.3e1_dp &
979 0 : /0.2e1_dp*rs_brhob*t464 + t150*t66*rs_brhob*t468
980 0 : t472 = 0.1e1_dp/t155
981 : epsilon_c_unif_brhob = -0.2e1_dp*t239*rs_brhob*t156 + t458* &
982 0 : t471*t472
983 0 : s_b_2rhob = 0.2e1_dp*s_b*s_brhob
984 0 : s_avg_2rhob = s_b_2rhob/0.2e1_dp
985 0 : e_lsda_c_brhob = epsilon_c_unif_brhob*my_rhob + epsilon_c_unif_b
986 0 : t480 = t339*s_avg_2rhob
987 0 : u_c_abrhob = gamma_c_ab*s_avg_2rhob*t162 - t337*t480
988 0 : t484 = t344*s_b_2
989 0 : t485 = t167**2
990 0 : t486 = 0.1e1_dp/t485
991 0 : 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 0 : e_lsda_c_brhob
994 0 : gc_abrhob = u_c_abrhob*t170 + u_c_ab*u_c_abrhob*c_cab_2
995 0 : gc_brhob = u_c_brhob*t176 + u_c_b*u_c_brhob*c_css_2
996 :
997 0 : 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 0 : gc_abrhob + e_lsda_c_brhob*gc_b + e_lsda_c_b*gc_brhob)
1001 0 : e_rb(ii) = e_rb(ii) + exc_rhob
1002 : END IF
1003 :
1004 0 : s_anorm_drhoa = t12
1005 : u_x_anorm_drhoa = 0.2e1_dp*t191*t16*s_anorm_drhoa - 0.2e1_dp &
1006 0 : *t196*t198*s_anorm_drhoa
1007 0 : gx_anorm_drhoa = u_x_anorm_drhoa*t18 + u_x_a*u_x_anorm_drhoa*c_x_2
1008 0 : s_a_2norm_drhoa = 0.2e1_dp*s_a*s_anorm_drhoa
1009 0 : s_avg_2norm_drhoa = s_a_2norm_drhoa/0.2e1_dp
1010 0 : t512 = t339*s_avg_2norm_drhoa
1011 0 : u_c_abnorm_drhoa = gamma_c_ab*s_avg_2norm_drhoa*t162 - t337*t512
1012 0 : t516 = t347*s_a_2norm_drhoa
1013 0 : 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 0 : u_c_abnorm_drhoa*c_cab_2
1016 : gc_anorm_drhoa = u_c_anorm_drhoa*t173 + u_c_a*u_c_anorm_drhoa &
1017 0 : *c_css_2
1018 :
1019 0 : 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 0 : (e_lsda_c_ab*gc_abnorm_drhoa + e_lsda_c_a*gc_anorm_drhoa)
1022 0 : e_ndra(ii) = e_ndra(ii) + exc_norm_drhoa
1023 : END IF
1024 :
1025 0 : s_bnorm_drhob = t25
1026 : u_x_bnorm_drhob = 0.2e1_dp*t370*t29*s_bnorm_drhob - 0.2e1_dp &
1027 0 : *t374*t376*s_bnorm_drhob
1028 0 : gx_bnorm_drhob = u_x_bnorm_drhob*t31 + u_x_b*u_x_bnorm_drhob*c_x_2
1029 0 : s_b_2norm_drhob = 0.2e1_dp*s_b*s_bnorm_drhob
1030 0 : s_avg_2norm_drhob = s_b_2norm_drhob/0.2e1_dp
1031 0 : t539 = t339*s_avg_2norm_drhob
1032 0 : u_c_abnorm_drhob = gamma_c_ab*s_avg_2norm_drhob*t162 - t337*t539
1033 0 : t543 = t486*s_b_2norm_drhob
1034 0 : 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 0 : u_c_abnorm_drhob*c_cab_2
1037 : gc_bnorm_drhob = u_c_bnorm_drhob*t176 + u_c_b*u_c_bnorm_drhob &
1038 0 : *c_css_2
1039 :
1040 0 : 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 0 : (e_lsda_c_ab*gc_abnorm_drhob + e_lsda_c_b*gc_bnorm_drhob)
1043 0 : e_ndrb(ii) = e_ndrb(ii) + exc_norm_drhob
1044 : END IF
1045 :
1046 0 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
1047 0 : t555 = t7**2
1048 0 : t560 = t186*my_rhoa
1049 0 : s_arhoarhoa = 0.28e2_dp/0.9e1_dp*my_norm_drhoa/t7/t560
1050 0 : t564 = s_arhoa**2
1051 0 : t568 = t194*s_a_2
1052 0 : t575 = t194*gamma_x
1053 0 : t576 = s_a_2**2
1054 0 : t577 = t575*t576
1055 0 : 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 0 : t577*t579*t564 - 0.2e1_dp*t196*t198*s_arhoarhoa
1059 0 : u_x_a1rhoa = u_x_arhoa
1060 0 : t600 = 0.1e1_dp/t207/rho
1061 0 : t601 = t33*t600
1062 0 : chirhoarhoa = -0.2e1_dp*t208 + 0.2e1_dp*t601
1063 0 : t605 = 0.3141592654e1_dp**2
1064 0 : t606 = 0.1e1_dp/t605
1065 0 : t608 = t207**2
1066 : rsrhoarhoa = -t4/t210/t36*t606/t608/0.18e2_dp + &
1067 0 : t4*t212*t600/0.6e1_dp
1068 0 : t619 = alpha_1_1*rsrhoa
1069 0 : t621 = t221*t235*t236
1070 0 : t626 = t40/t220/t51
1071 0 : t627 = t235**2
1072 0 : t631 = 0.1e1_dp/t46
1073 0 : t632 = beta_1_1*t631
1074 0 : t633 = rsrhoa**2
1075 0 : t639 = beta_3_1*t223
1076 0 : t644 = t48**2
1077 0 : t646 = rs**2
1078 0 : t647 = 0.1e1_dp/t646
1079 0 : t659 = t220**2
1080 0 : t661 = t40/t659
1081 0 : t662 = t55**2
1082 0 : 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 : 0.2e1_dp
1090 0 : e_c_u_01rhoa = e_c_u_0rhoa
1091 0 : t671 = alpha_1_2*rsrhoa
1092 0 : t673 = t244*t256*t257
1093 0 : t678 = t60/t243/t69
1094 0 : t679 = t256**2
1095 0 : t683 = beta_1_2*t631
1096 0 : t689 = beta_3_2*t223
1097 0 : t694 = t66**2
1098 0 : t707 = t243**2
1099 0 : t709 = t60/t707
1100 0 : t710 = t73**2
1101 0 : t711 = 0.1e1_dp/t710
1102 0 : t719 = alpha_1_3*rsrhoa
1103 0 : t721 = t265*t277*t278
1104 0 : t726 = t78/t264/t87
1105 0 : t727 = t277**2
1106 0 : t731 = beta_1_3*t631
1107 0 : t737 = beta_3_3*t223
1108 0 : t742 = t84**2
1109 0 : t755 = t264**2
1110 0 : t757 = t78/t755
1111 0 : t758 = t91**2
1112 0 : t759 = 0.1e1_dp/t758
1113 0 : alpha_c1rhoa = alpha_crhoa
1114 0 : t764 = t99**2
1115 0 : t765 = 0.1e1_dp/t764
1116 0 : t766 = chirhoa**2
1117 0 : t771 = t102**2
1118 0 : 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 : 0.4e1_dp/0.3e1_dp*t102*chirhoarhoa)*t97
1122 0 : f1rhoa = frhoa
1123 0 : t790 = alpha_c1rhoa*f
1124 0 : t793 = alpha_c*f1rhoa
1125 0 : t796 = t106*t107
1126 0 : t811 = e_c_u_1rhoa - e_c_u_01rhoa
1127 0 : t818 = t811*f
1128 0 : 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 0 : t766 + 0.4e1_dp*t113*t289*chirhoarhoa
1140 : epsilon_c_unif1rhoa = e_c_u_01rhoa + t790*t110 + t793*t110 - &
1141 0 : t293 + t818*t108 + t821*t108 + t301
1142 0 : t838 = t186**2
1143 : rs_arhoarhoa = -t4/t302/t116*t606/t838/0.18e2_dp + &
1144 0 : t4*t304/t560/0.6e1_dp
1145 0 : t858 = t327**2
1146 0 : t864 = rs_arhoa**2
1147 0 : t876 = rs_a**2
1148 0 : t877 = 0.1e1_dp/t876
1149 0 : t889 = t312**2
1150 0 : t892 = t133**2
1151 0 : epsilon_c_unif_a1rhoa = epsilon_c_unif_arhoa
1152 0 : s_a_2rhoarhoa = 0.2e1_dp*t564 + 0.2e1_dp*s_a*s_arhoarhoa
1153 0 : s_a_21rhoa = s_a_2rhoa
1154 0 : s_avg_2rhoarhoa = s_a_2rhoarhoa/0.2e1_dp
1155 0 : 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 0 : + epsilon_c_unif_a1rhoa
1165 0 : e_lsda_c_a1rhoa = epsilon_c_unif_a1rhoa*my_rhoa + epsilon_c_unif_a
1166 0 : t906 = t336*s_avg_2rhoa
1167 0 : t907 = t339*s_avg_21rhoa
1168 0 : t911 = t336*gamma_c_ab*s_avg_2
1169 0 : t913 = 0.1e1_dp/t338/t161
1170 0 : 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 0 : s_avg_2rhoarhoa
1174 0 : u_c_ab1rhoa = gamma_c_ab*s_avg_21rhoa*t162 - t337*t907
1175 0 : t925 = t344*s_a_2rhoa
1176 0 : t926 = t347*s_a_21rhoa
1177 0 : t929 = t344*gamma_c_ss
1178 0 : t930 = t929*s_a_2
1179 0 : t932 = 0.1e1_dp/t346/t164
1180 0 : 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 0 : s_a_2rhoarhoa
1184 0 : u_c_a1rhoa = gamma_c_ss*s_a_21rhoa*t165 - t345*t926
1185 0 : 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 0 : u_c_a1rhoa*c_css_2 + u_c_a*u_c_arhoarhoa*c_css_2))
1210 0 : e_ra_ra(ii) = e_ra_ra(ii) + exc_rhoa_rhoa
1211 : END IF
1212 :
1213 0 : chirhoarhob = 0.2e1_dp*t601
1214 0 : rsrhoarhob = rsrhoarhoa
1215 0 : t974 = t221*t396*t236
1216 0 : t976 = alpha_1_1*rsrhob
1217 0 : t981 = rsrhoa*rsrhob
1218 0 : 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 : 0.2e1_dp
1226 0 : t1012 = t244*t410*t257
1227 0 : t1014 = alpha_1_2*rsrhob
1228 0 : t1047 = t265*t424*t278
1229 0 : 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 0 : t97
1234 0 : 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 : 0.4e1_dp*t113*t289*chirhoarhob
1246 : u_c_abrhoarhob = -0.2e1_dp*t906*t480 + 0.2e1_dp*t911*t914 &
1247 0 : *s_avg_2rhob
1248 :
1249 0 : 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 0 : u_c_abrhob*c_cab_2 + u_c_ab*u_c_abrhoarhob*c_cab_2))
1264 0 : e_ra_rb(ii) = e_ra_rb(ii) + exc_rhoa_rhob
1265 : END IF
1266 :
1267 0 : t1152 = t20**2
1268 0 : t1157 = t365*my_rhob
1269 0 : s_brhobrhob = 0.28e2_dp/0.9e1_dp*my_norm_drhob/t20/t1157
1270 0 : t1161 = s_brhob**2
1271 0 : t1165 = t194*s_b_2
1272 0 : t1172 = s_b_2**2
1273 0 : t1173 = t575*t1172
1274 0 : 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 0 : s_brhobrhob
1279 0 : u_x_b1rhob = u_x_brhob
1280 0 : chirhobrhob = 0.2e1_dp*t208 + 0.2e1_dp*t601
1281 0 : rsrhobrhob = rsrhoarhob
1282 0 : t1201 = t396**2
1283 0 : 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 0 : t1201*t663*t42/0.2e1_dp
1291 0 : e_c_u_01rhob = e_c_u_0rhob
1292 0 : t1236 = t410**2
1293 0 : t1270 = t424**2
1294 0 : alpha_c1rhob = alpha_crhob
1295 0 : 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 : 0.4e1_dp/0.3e1_dp*t102*chirhobrhob)*t97
1299 0 : f1rhob = frhob
1300 0 : t1321 = alpha_c1rhob*f
1301 0 : t1324 = alpha_c*f1rhob
1302 0 : t1341 = e_c_u_1rhob - e_c_u_01rhob
1303 0 : t1348 = t1341*f
1304 0 : 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 0 : *t107*t1299 + 0.4e1_dp*t113*t289*chirhobrhob
1316 : epsilon_c_unif1rhob = e_c_u_01rhob + t1321*t110 + t1324*t110 - &
1317 0 : t437 + t1348*t108 + t1351*t108 + t445
1318 0 : t1368 = t365**2
1319 : rs_brhobrhob = -t4/t446/t138*t606/t1368/0.18e2_dp &
1320 0 : + t4*t448/t1157/0.6e1_dp
1321 0 : t1388 = t471**2
1322 0 : t1394 = rs_brhob**2
1323 0 : t1406 = rs_b**2
1324 0 : t1407 = 0.1e1_dp/t1406
1325 0 : t1419 = t456**2
1326 0 : t1422 = t155**2
1327 0 : epsilon_c_unif_b1rhob = epsilon_c_unif_brhob
1328 0 : s_b_2rhobrhob = 0.2e1_dp*t1161 + 0.2e1_dp*s_b*s_brhobrhob
1329 0 : s_b_21rhob = s_b_2rhob
1330 0 : s_avg_2rhobrhob = s_b_2rhobrhob/0.2e1_dp
1331 0 : 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 0 : my_rhob + epsilon_c_unif_brhob + epsilon_c_unif_b1rhob
1341 0 : e_lsda_c_b1rhob = epsilon_c_unif_b1rhob*my_rhob + epsilon_c_unif_b
1342 0 : t1436 = t336*s_avg_2rhob
1343 0 : t1437 = t339*s_avg_21rhob
1344 0 : 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 0 : *s_avg_2rhobrhob
1348 0 : u_c_ab1rhob = gamma_c_ab*s_avg_21rhob*t162 - t337*t1437
1349 0 : t1451 = t344*s_b_2rhob
1350 0 : t1452 = t486*s_b_21rhob
1351 0 : t1455 = t929*s_b_2
1352 0 : t1457 = 0.1e1_dp/t485/t167
1353 0 : 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 0 : *s_b_2rhobrhob
1357 0 : u_c_b1rhob = gamma_c_ss*s_b_21rhob*t168 - t484*t1452
1358 :
1359 0 : 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 0 : *c_css_2))
1385 0 : e_rb_rb(ii) = e_rb_rb(ii) + exc_rhob_rhob
1386 : END IF
1387 :
1388 0 : 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 : - 0.2e1_dp*t196*t198*s_arhoanorm_drhoa
1393 : s_a_2rhoanorm_drhoa = 0.2e1_dp*s_anorm_drhoa*s_arhoa + &
1394 0 : 0.2e1_dp*s_a*s_arhoanorm_drhoa
1395 0 : 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 0 : - 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 0 : t345*t347*s_a_2rhoanorm_drhoa
1402 :
1403 0 : 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 0 : *c_css_2 + u_c_a*u_c_arhoanorm_drhoa*c_css_2))
1413 0 : 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 0 : *t1440*s_avg_2norm_drhoa
1418 :
1419 0 : 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 0 : u_c_abrhobnorm_drhoa*c_cab_2))
1424 0 : e_rb_ndra(ii) = e_rb_ndra(ii) + exc_rhob_norm_drhoa
1425 : END IF
1426 :
1427 0 : t1571 = s_anorm_drhoa**2
1428 : u_x_anorm_drhoanorm_drhoa = 0.2e1_dp*gamma_x*t1571*t16 - &
1429 0 : 0.10e2_dp*t568*t198*t1571 + 0.8e1_dp*t577*t579*t1571
1430 0 : s_a_2norm_drhoanorm_drhoa = 0.2e1_dp*t1571
1431 0 : s_a_21norm_drhoa = s_a_2norm_drhoa
1432 0 : s_avg_2norm_drhoanorm_drhoa = s_a_2norm_drhoanorm_drhoa/0.2e1_dp
1433 0 : s_avg_21norm_drhoa = s_a_21norm_drhoa/0.2e1_dp
1434 0 : t1589 = t336*s_avg_2norm_drhoa
1435 0 : t1590 = t339*s_avg_21norm_drhoa
1436 0 : 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 0 : s_avg_2norm_drhoanorm_drhoa
1441 0 : 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 0 : s_a_2norm_drhoanorm_drhoa
1446 :
1447 0 : 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 0 : t1605)*c_css_2 + u_c_a*u_c_anorm_drhoanorm_drhoa*c_css_2))
1457 0 : 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 0 : t914*s_avg_2norm_drhob
1462 :
1463 0 : 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 0 : u_c_abrhoanorm_drhob*c_cab_2))
1468 0 : e_ra_ndrb(ii) = e_ra_ndrb(ii) + exc_rhoa_norm_drhob
1469 : END IF
1470 :
1471 0 : 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 0 : 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 : 0.2e1_dp*s_b*s_brhobnorm_drhob
1478 0 : 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 0 : 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 0 : - t484*t486*s_b_2rhobnorm_drhob
1485 :
1486 0 : 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 0 : *c_css_2 + u_c_b*u_c_brhobnorm_drhob*c_css_2))
1496 0 : 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 0 : t911*t1593*s_avg_2norm_drhob
1501 :
1502 0 : 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 0 : c_cab_2)
1507 0 : e_ndra_ndrb(ii) = e_ndra_ndrb(ii) + exc_norm_drhoa_norm_drhob
1508 : END IF
1509 :
1510 0 : t1719 = s_bnorm_drhob**2
1511 : u_x_bnorm_drhobnorm_drhob = 0.2e1_dp*gamma_x*t1719*t29 - &
1512 0 : 0.10e2_dp*t1165*t376*t1719 + 0.8e1_dp*t1173*t1175*t1719
1513 0 : s_b_2norm_drhobnorm_drhob = 0.2e1_dp*t1719
1514 0 : s_b_21norm_drhob = s_b_2norm_drhob
1515 0 : s_avg_2norm_drhobnorm_drhob = s_b_2norm_drhobnorm_drhob/0.2e1_dp
1516 0 : s_avg_21norm_drhob = s_b_21norm_drhob/0.2e1_dp
1517 0 : 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 0 : s_avg_2norm_drhobnorm_drhob
1523 0 : 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 0 : s_b_2norm_drhobnorm_drhob
1528 :
1529 0 : 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 0 : t1753)*c_css_2 + u_c_b*u_c_bnorm_drhobnorm_drhob*c_css_2))
1539 0 : 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 0 : 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 607 : 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 607 : p_1 = 0.10e1_dp
1648 607 : A_1 = 0.31091e-1_dp
1649 607 : alpha_1_1 = 0.21370e0_dp
1650 607 : beta_1_1 = 0.75957e1_dp
1651 607 : beta_2_1 = 0.35876e1_dp
1652 607 : beta_3_1 = 0.16382e1_dp
1653 607 : beta_4_1 = 0.49294e0_dp
1654 607 : p_2 = 0.10e1_dp
1655 607 : A_2 = 0.15545e-1_dp
1656 607 : alpha_1_2 = 0.20548e0_dp
1657 607 : beta_1_2 = 0.141189e2_dp
1658 607 : beta_2_2 = 0.61977e1_dp
1659 607 : beta_3_2 = 0.33662e1_dp
1660 607 : beta_4_2 = 0.62517e0_dp
1661 607 : p_3 = 0.10e1_dp
1662 607 : A_3 = 0.16887e-1_dp
1663 607 : alpha_1_3 = 0.11125e0_dp
1664 607 : beta_1_3 = 0.10357e2_dp
1665 607 : beta_2_3 = 0.36231e1_dp
1666 607 : beta_3_3 = 0.88026e0_dp
1667 607 : beta_4_3 = 0.49671e0_dp
1668 :
1669 607 : coeffs = b97_coeffs(param)
1670 607 : c_x_0 = coeffs(1)
1671 607 : c_x_1 = coeffs(2)
1672 607 : c_x_2 = coeffs(3)
1673 607 : c_cab_0 = coeffs(4)
1674 607 : c_cab_1 = coeffs(5)
1675 607 : c_cab_2 = coeffs(6)
1676 607 : c_css_0 = coeffs(7)
1677 607 : c_css_1 = coeffs(8)
1678 607 : c_css_2 = coeffs(9)
1679 :
1680 607 : scale_c = scale_c_in
1681 607 : scale_x = scale_x_in
1682 607 : IF (scale_x == -1.0_dp) scale_x = coeffs(10)
1683 :
1684 607 : gamma_x = 0.004_dp
1685 607 : gamma_c_ss = 0.2_dp
1686 607 : gamma_c_ab = 0.006_dp
1687 :
1688 607 : SELECT CASE (grad_deriv)
1689 : CASE default
1690 607 : t1 = 3**(0.1e1_dp/0.3e1_dp)
1691 607 : t2 = 4**(0.1e1_dp/0.3e1_dp)
1692 607 : t3 = t2**2
1693 607 : t4 = t1*t3
1694 607 : t6 = (0.1e1_dp/pi)**(0.1e1_dp/0.3e1_dp)
1695 607 : !$OMP DO
1696 : DO ii = 1, npoints
1697 15268288 : my_rhoa = 0.5_dp*MAX(rho_tot(ii), 0.0_dp)
1698 15268288 : my_rhob = my_rhoa
1699 15268288 : rho = my_rhoa + my_rhob
1700 15268288 : IF (rho > epsilon_rho) THEN
1701 8586461 : my_rhoa = MAX(my_rhoa, 0.5_dp*epsilon_rho)
1702 8586461 : my_rhob = MAX(my_rhob, 0.5_dp*epsilon_rho)
1703 8586461 : rho = my_rhoa + my_rhob
1704 8586461 : my_norm_drhoa = 0.5_dp*norm_drho(ii)
1705 8586461 : my_norm_drhob = 0.5_dp*norm_drho(ii)
1706 8586461 : t7 = my_rhoa**(0.1e1_dp/0.3e1_dp)
1707 8586461 : t8 = t7*my_rhoa
1708 8586461 : e_lsda_x_a = -0.3e1_dp/0.8e1_dp*t4*t6*t8
1709 8586461 : t12 = 0.1e1_dp/t8
1710 8586461 : s_a = my_norm_drhoa*t12
1711 8586461 : t13 = s_a**2
1712 8586461 : t14 = gamma_x*t13
1713 8586461 : t15 = 0.1e1_dp + t14
1714 8586461 : t16 = 0.1e1_dp/t15
1715 8586461 : u_x_a = t14*t16
1716 8586461 : t18 = c_x_1 + u_x_a*c_x_2
1717 8586461 : gx_a = c_x_0 + u_x_a*t18
1718 8586461 : t20 = my_rhob**(0.1e1_dp/0.3e1_dp)
1719 8586461 : t21 = t20*my_rhob
1720 8586461 : e_lsda_x_b = -0.3e1_dp/0.8e1_dp*t4*t6*t21
1721 8586461 : t25 = 0.1e1_dp/t21
1722 8586461 : s_b = my_norm_drhob*t25
1723 8586461 : t26 = s_b**2
1724 8586461 : t27 = gamma_x*t26
1725 8586461 : t28 = 0.1e1_dp + t27
1726 8586461 : t29 = 0.1e1_dp/t28
1727 8586461 : u_x_b = t27*t29
1728 8586461 : t31 = c_x_1 + u_x_b*c_x_2
1729 8586461 : gx_b = c_x_0 + u_x_b*t31
1730 8586461 : t33 = my_rhoa - my_rhob
1731 8586461 : t34 = 0.1e1_dp/rho
1732 8586461 : chi = t33*t34
1733 8586461 : t35 = 0.1e1_dp/pi
1734 8586461 : t36 = t35*t34
1735 8586461 : t37 = t36**(0.1e1_dp/0.3e1_dp)
1736 8586461 : rs = t4*t37/0.4e1_dp
1737 8586461 : t40 = 0.1e1_dp + alpha_1_1*rs
1738 8586461 : t42 = 0.1e1_dp/A_1
1739 8586461 : t43 = SQRT(rs)
1740 8586461 : t46 = t43*rs
1741 8586461 : t48 = p_1 + 0.1e1_dp
1742 8586461 : t49 = rs**t48
1743 8586461 : t50 = beta_4_1*t49
1744 8586461 : t51 = beta_1_1*t43 + beta_2_1*rs + beta_3_1*t46 + t50
1745 8586461 : t55 = 0.1e1_dp + t42/t51/0.2e1_dp
1746 8586461 : t56 = LOG(t55)
1747 8586461 : e_c_u_0 = -0.2e1_dp*A_1*t40*t56
1748 8586461 : t60 = 0.1e1_dp + alpha_1_2*rs
1749 8586461 : t62 = 0.1e1_dp/A_2
1750 8586461 : t66 = p_2 + 0.1e1_dp
1751 8586461 : t67 = rs**t66
1752 8586461 : t68 = beta_4_2*t67
1753 8586461 : t69 = beta_1_2*t43 + beta_2_2*rs + beta_3_2*t46 + t68
1754 8586461 : t73 = 0.1e1_dp + t62/t69/0.2e1_dp
1755 8586461 : t74 = LOG(t73)
1756 8586461 : t78 = 0.1e1_dp + alpha_1_3*rs
1757 8586461 : t80 = 0.1e1_dp/A_3
1758 8586461 : t84 = p_3 + 0.1e1_dp
1759 8586461 : t85 = rs**t84
1760 8586461 : t86 = beta_4_3*t85
1761 8586461 : t87 = beta_1_3*t43 + beta_2_3*rs + beta_3_3*t46 + t86
1762 8586461 : t91 = 0.1e1_dp + t80/t87/0.2e1_dp
1763 8586461 : t92 = LOG(t91)
1764 8586461 : alpha_c = 0.2e1_dp*A_3*t78*t92
1765 8586461 : t94 = 2**(0.1e1_dp/0.3e1_dp)
1766 8586461 : t97 = 1/(2*t94 - 2)
1767 8586461 : t98 = 0.1e1_dp + chi
1768 8586461 : t99 = t98**(0.1e1_dp/0.3e1_dp)
1769 8586461 : t101 = 0.1e1_dp - chi
1770 8586461 : t102 = t101**(0.1e1_dp/0.3e1_dp)
1771 8586461 : f = (t99*t98 + t102*t101 - 0.2e1_dp)*t97
1772 8586461 : t105 = alpha_c*f
1773 8586461 : t106 = 0.9e1_dp/0.8e1_dp/t97
1774 8586461 : t107 = chi**2
1775 8586461 : t108 = t107**2
1776 8586461 : t110 = t106*(0.1e1_dp - t108)
1777 8586461 : t112 = -0.2e1_dp*A_2*t60*t74 - e_c_u_0
1778 8586461 : t113 = t112*f
1779 8586461 : epsilon_c_unif = e_c_u_0 + t105*t110 + t113*t108
1780 8586461 : t116 = t35/my_rhoa
1781 8586461 : t117 = t116**(0.1e1_dp/0.3e1_dp)
1782 8586461 : rs_a = t4*t117/0.4e1_dp
1783 8586461 : t120 = 0.1e1_dp + alpha_1_2*rs_a
1784 8586461 : t122 = SQRT(rs_a)
1785 8586461 : t125 = t122*rs_a
1786 8586461 : t127 = rs_a**t66
1787 8586461 : t128 = beta_4_2*t127
1788 8586461 : t129 = beta_1_2*t122 + beta_2_2*rs_a + beta_3_2*t125 + t128
1789 8586461 : t133 = 0.1e1_dp + t62/t129/0.2e1_dp
1790 8586461 : t134 = LOG(t133)
1791 8586461 : epsilon_c_unif_a = -0.2e1_dp*A_2*t120*t134
1792 8586461 : t138 = t35/my_rhob
1793 8586461 : t139 = t138**(0.1e1_dp/0.3e1_dp)
1794 8586461 : rs_b = t4*t139/0.4e1_dp
1795 8586461 : t142 = 0.1e1_dp + alpha_1_2*rs_b
1796 8586461 : t144 = SQRT(rs_b)
1797 8586461 : t147 = t144*rs_b
1798 8586461 : t149 = rs_b**t66
1799 8586461 : t150 = beta_4_2*t149
1800 8586461 : t151 = beta_1_2*t144 + beta_2_2*rs_b + beta_3_2*t147 + t150
1801 8586461 : t155 = 0.1e1_dp + t62/t151/0.2e1_dp
1802 8586461 : t156 = LOG(t155)
1803 8586461 : epsilon_c_unif_b = -0.2e1_dp*A_2*t142*t156
1804 8586461 : s_a_2 = t13
1805 8586461 : s_b_2 = t26
1806 8586461 : s_avg_2 = s_a_2/0.2e1_dp + s_b_2/0.2e1_dp
1807 8586461 : e_lsda_c_a = epsilon_c_unif_a*my_rhoa
1808 8586461 : e_lsda_c_b = epsilon_c_unif_b*my_rhob
1809 8586461 : t160 = gamma_c_ab*s_avg_2
1810 8586461 : t161 = 0.1e1_dp + t160
1811 8586461 : t162 = 0.1e1_dp/t161
1812 8586461 : u_c_ab = t160*t162
1813 8586461 : t163 = gamma_c_ss*s_a_2
1814 8586461 : t164 = 0.1e1_dp + t163
1815 8586461 : t165 = 0.1e1_dp/t164
1816 8586461 : u_c_a = t163*t165
1817 8586461 : t166 = gamma_c_ss*s_b_2
1818 8586461 : t167 = 0.1e1_dp + t166
1819 8586461 : t168 = 0.1e1_dp/t167
1820 8586461 : u_c_b = t166*t168
1821 8586461 : e_lsda_c_ab = epsilon_c_unif*rho - e_lsda_c_a - e_lsda_c_b
1822 8586461 : t170 = c_cab_1 + u_c_ab*c_cab_2
1823 8586461 : gc_ab = c_cab_0 + u_c_ab*t170
1824 8586461 : t173 = c_css_1 + u_c_a*c_css_2
1825 8586461 : gc_a = c_css_0 + u_c_a*t173
1826 8586461 : t176 = c_css_1 + u_c_b*c_css_2
1827 8586461 : gc_b = c_css_0 + u_c_b*t176
1828 :
1829 8586461 : IF (grad_deriv >= 0) THEN
1830 : exc = scale_x*(e_lsda_x_a*gx_a + e_lsda_x_b*gx_b) + scale_c &
1831 8586461 : *(e_lsda_c_ab*gc_ab + e_lsda_c_a*gc_a + e_lsda_c_b*gc_b)
1832 8586461 : e_0(ii) = e_0(ii) + exc
1833 : END IF
1834 :
1835 8586461 : IF (grad_deriv /= 0) THEN
1836 :
1837 8586461 : e_lsda_x_arhoa = -t4*t6*t7/0.2e1_dp
1838 8586461 : t186 = my_rhoa**2
1839 8586461 : t188 = 0.1e1_dp/t7/t186
1840 8586461 : s_arhoa = -0.4e1_dp/0.3e1_dp*my_norm_drhoa*t188
1841 8586461 : t191 = gamma_x*s_a
1842 8586461 : t192 = t16*s_arhoa
1843 8586461 : t194 = gamma_x**2
1844 8586461 : t196 = t194*s_a_2*s_a
1845 8586461 : t197 = t15**2
1846 8586461 : t198 = 0.1e1_dp/t197
1847 8586461 : t199 = t198*s_arhoa
1848 8586461 : u_x_arhoa = 0.2e1_dp*t191*t192 - 0.2e1_dp*t196*t199
1849 8586461 : gx_arhoa = u_x_arhoa*t18 + u_x_a*u_x_arhoa*c_x_2
1850 8586461 : t207 = rho**2
1851 8586461 : t208 = 0.1e1_dp/t207
1852 8586461 : t209 = t33*t208
1853 8586461 : chirhoa = t34 - t209
1854 8586461 : t210 = t37**2
1855 8586461 : t212 = 0.1e1_dp/t210*t35
1856 8586461 : rsrhoa = -t4*t212*t208/0.12e2_dp
1857 8586461 : t216 = A_1*alpha_1_1
1858 8586461 : t220 = t51**2
1859 8586461 : t221 = 0.1e1_dp/t220
1860 8586461 : t222 = t40*t221
1861 8586461 : t223 = 0.1e1_dp/t43
1862 8586461 : t224 = beta_1_1*t223
1863 8586461 : t228 = beta_3_1*t43
1864 8586461 : t232 = 0.1e1_dp/rs
1865 : t235 = t224*rsrhoa/0.2e1_dp + beta_2_1*rsrhoa + 0.3e1_dp/ &
1866 8586461 : 0.2e1_dp*t228*rsrhoa + t50*t48*rsrhoa*t232
1867 8586461 : t236 = 0.1e1_dp/t55
1868 8586461 : t237 = t235*t236
1869 8586461 : e_c_u_0rhoa = -0.2e1_dp*t216*rsrhoa*t56 + t222*t237
1870 8586461 : t239 = A_2*alpha_1_2
1871 8586461 : t243 = t69**2
1872 8586461 : t244 = 0.1e1_dp/t243
1873 8586461 : t245 = t60*t244
1874 8586461 : t246 = beta_1_2*t223
1875 8586461 : t250 = beta_3_2*t43
1876 : t256 = t246*rsrhoa/0.2e1_dp + beta_2_2*rsrhoa + 0.3e1_dp/ &
1877 8586461 : 0.2e1_dp*t250*rsrhoa + t68*t66*rsrhoa*t232
1878 8586461 : t257 = 0.1e1_dp/t73
1879 8586461 : t258 = t256*t257
1880 8586461 : e_c_u_1rhoa = -0.2e1_dp*t239*rsrhoa*t74 + t245*t258
1881 8586461 : t260 = A_3*alpha_1_3
1882 8586461 : t264 = t87**2
1883 8586461 : t265 = 0.1e1_dp/t264
1884 8586461 : t266 = t78*t265
1885 8586461 : t267 = beta_1_3*t223
1886 8586461 : t271 = beta_3_3*t43
1887 : t277 = t267*rsrhoa/0.2e1_dp + beta_2_3*rsrhoa + 0.3e1_dp/ &
1888 8586461 : 0.2e1_dp*t271*rsrhoa + t86*t84*rsrhoa*t232
1889 8586461 : t278 = 0.1e1_dp/t91
1890 8586461 : t279 = t277*t278
1891 8586461 : 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 8586461 : *t102*chirhoa)*t97
1894 8586461 : t285 = alpha_crhoa*f
1895 8586461 : t287 = alpha_c*frhoa
1896 8586461 : t289 = t107*chi
1897 8586461 : t290 = t106*t289
1898 8586461 : t291 = t290*chirhoa
1899 8586461 : t293 = 0.4e1_dp*t105*t291
1900 8586461 : t294 = e_c_u_1rhoa - e_c_u_0rhoa
1901 8586461 : t295 = t294*f
1902 8586461 : t297 = t112*frhoa
1903 8586461 : t299 = t289*chirhoa
1904 8586461 : t301 = 0.4e1_dp*t113*t299
1905 : epsilon_c_unifrhoa = e_c_u_0rhoa + t285*t110 + t287*t110 - &
1906 8586461 : t293 + t295*t108 + t297*t108 + t301
1907 8586461 : t302 = t117**2
1908 8586461 : t304 = 0.1e1_dp/t302*t35
1909 8586461 : rs_arhoa = -t4*t304/t186/0.12e2_dp
1910 8586461 : t312 = t129**2
1911 8586461 : t313 = 0.1e1_dp/t312
1912 8586461 : t314 = t120*t313
1913 8586461 : t315 = 0.1e1_dp/t122
1914 8586461 : t316 = beta_1_2*t315
1915 8586461 : t320 = beta_3_2*t122
1916 8586461 : t324 = 0.1e1_dp/rs_a
1917 : t327 = t316*rs_arhoa/0.2e1_dp + beta_2_2*rs_arhoa + 0.3e1_dp &
1918 8586461 : /0.2e1_dp*t320*rs_arhoa + t128*t66*rs_arhoa*t324
1919 8586461 : t328 = 0.1e1_dp/t133
1920 : epsilon_c_unif_arhoa = -0.2e1_dp*t239*rs_arhoa*t134 + t314* &
1921 8586461 : t327*t328
1922 8586461 : s_a_2rhoa = 0.2e1_dp*s_a*s_arhoa
1923 8586461 : s_avg_2rhoa = s_a_2rhoa/0.2e1_dp
1924 8586461 : e_lsda_c_arhoa = epsilon_c_unif_arhoa*my_rhoa + epsilon_c_unif_a
1925 8586461 : t336 = gamma_c_ab**2
1926 8586461 : t337 = t336*s_avg_2
1927 8586461 : t338 = t161**2
1928 8586461 : t339 = 0.1e1_dp/t338
1929 8586461 : u_c_abrhoa = gamma_c_ab*s_avg_2rhoa*t162 - t337*t339*s_avg_2rhoa
1930 8586461 : t344 = gamma_c_ss**2
1931 8586461 : t345 = t344*s_a_2
1932 8586461 : t346 = t164**2
1933 8586461 : t347 = 0.1e1_dp/t346
1934 8586461 : 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 8586461 : e_lsda_c_arhoa
1937 8586461 : gc_abrhoa = u_c_abrhoa*t170 + u_c_ab*u_c_abrhoa*c_cab_2
1938 8586461 : gc_arhoa = u_c_arhoa*t173 + u_c_a*u_c_arhoa*c_css_2
1939 :
1940 8586461 : 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 8586461 : gc_abrhoa + e_lsda_c_arhoa*gc_a + e_lsda_c_a*gc_arhoa)
1944 8586461 : e_r(ii) = e_r(ii) + 0.5_dp*exc_rhoa
1945 : END IF
1946 :
1947 8586461 : e_lsda_x_brhob = -t4*t6*t20/0.2e1_dp
1948 8586461 : t365 = my_rhob**2
1949 8586461 : t367 = 0.1e1_dp/t20/t365
1950 8586461 : s_brhob = -0.4e1_dp/0.3e1_dp*my_norm_drhob*t367
1951 8586461 : t370 = gamma_x*s_b
1952 8586461 : t371 = t29*s_brhob
1953 8586461 : t374 = t194*s_b_2*s_b
1954 8586461 : t375 = t28**2
1955 8586461 : t376 = 0.1e1_dp/t375
1956 8586461 : t377 = t376*s_brhob
1957 8586461 : u_x_brhob = 0.2e1_dp*t370*t371 - 0.2e1_dp*t374*t377
1958 8586461 : gx_brhob = u_x_brhob*t31 + u_x_b*u_x_brhob*c_x_2
1959 8586461 : chirhob = -t34 - t209
1960 8586461 : rsrhob = rsrhoa
1961 : t396 = t224*rsrhob/0.2e1_dp + beta_2_1*rsrhob + 0.3e1_dp/ &
1962 8586461 : 0.2e1_dp*t228*rsrhob + t50*t48*rsrhob*t232
1963 8586461 : 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 8586461 : 0.2e1_dp*t250*rsrhob + t68*t66*rsrhob*t232
1966 8586461 : 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 8586461 : 0.2e1_dp*t271*rsrhob + t86*t84*rsrhob*t232
1969 8586461 : 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 8586461 : *t102*chirhob)*t97
1972 8586461 : t431 = alpha_crhob*f
1973 8586461 : t433 = alpha_c*frhob
1974 8586461 : t435 = t290*chirhob
1975 8586461 : t437 = 0.4e1_dp*t105*t435
1976 8586461 : t438 = e_c_u_1rhob - e_c_u_0rhob
1977 8586461 : t439 = t438*f
1978 8586461 : t441 = t112*frhob
1979 8586461 : t443 = t289*chirhob
1980 8586461 : t445 = 0.4e1_dp*t113*t443
1981 : epsilon_c_unifrhob = e_c_u_0rhob + t431*t110 + t433*t110 - &
1982 8586461 : t437 + t439*t108 + t441*t108 + t445
1983 8586461 : t446 = t139**2
1984 8586461 : t448 = 0.1e1_dp/t446*t35
1985 8586461 : rs_brhob = -t4*t448/t365/0.12e2_dp
1986 8586461 : t456 = t151**2
1987 8586461 : t457 = 0.1e1_dp/t456
1988 8586461 : t458 = t142*t457
1989 8586461 : t459 = 0.1e1_dp/t144
1990 8586461 : t460 = beta_1_2*t459
1991 8586461 : t464 = beta_3_2*t144
1992 8586461 : t468 = 0.1e1_dp/rs_b
1993 : t471 = t460*rs_brhob/0.2e1_dp + beta_2_2*rs_brhob + 0.3e1_dp &
1994 8586461 : /0.2e1_dp*rs_brhob*t464 + t150*t66*rs_brhob*t468
1995 8586461 : t472 = 0.1e1_dp/t155
1996 : epsilon_c_unif_brhob = -0.2e1_dp*t239*rs_brhob*t156 + t458* &
1997 8586461 : t471*t472
1998 8586461 : s_b_2rhob = 0.2e1_dp*s_b*s_brhob
1999 8586461 : s_avg_2rhob = s_b_2rhob/0.2e1_dp
2000 8586461 : e_lsda_c_brhob = epsilon_c_unif_brhob*my_rhob + epsilon_c_unif_b
2001 8586461 : t480 = t339*s_avg_2rhob
2002 8586461 : u_c_abrhob = gamma_c_ab*s_avg_2rhob*t162 - t337*t480
2003 8586461 : t484 = t344*s_b_2
2004 8586461 : t485 = t167**2
2005 8586461 : t486 = 0.1e1_dp/t485
2006 8586461 : 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 8586461 : e_lsda_c_brhob
2009 8586461 : gc_abrhob = u_c_abrhob*t170 + u_c_ab*u_c_abrhob*c_cab_2
2010 8586461 : gc_brhob = u_c_brhob*t176 + u_c_b*u_c_brhob*c_css_2
2011 :
2012 8586461 : 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 8586461 : gc_abrhob + e_lsda_c_brhob*gc_b + e_lsda_c_b*gc_brhob)
2016 8586461 : e_r(ii) = e_r(ii) + 0.5_dp*exc_rhob
2017 : END IF
2018 :
2019 8586461 : s_anorm_drhoa = t12
2020 : u_x_anorm_drhoa = 0.2e1_dp*t191*t16*s_anorm_drhoa - 0.2e1_dp &
2021 8586461 : *t196*t198*s_anorm_drhoa
2022 8586461 : gx_anorm_drhoa = u_x_anorm_drhoa*t18 + u_x_a*u_x_anorm_drhoa*c_x_2
2023 8586461 : s_a_2norm_drhoa = 0.2e1_dp*s_a*s_anorm_drhoa
2024 8586461 : s_avg_2norm_drhoa = s_a_2norm_drhoa/0.2e1_dp
2025 8586461 : t512 = t339*s_avg_2norm_drhoa
2026 8586461 : u_c_abnorm_drhoa = gamma_c_ab*s_avg_2norm_drhoa*t162 - t337*t512
2027 8586461 : t516 = t347*s_a_2norm_drhoa
2028 8586461 : 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 8586461 : u_c_abnorm_drhoa*c_cab_2
2031 : gc_anorm_drhoa = u_c_anorm_drhoa*t173 + u_c_a*u_c_anorm_drhoa &
2032 8586461 : *c_css_2
2033 :
2034 8586461 : 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 8586461 : (e_lsda_c_ab*gc_abnorm_drhoa + e_lsda_c_a*gc_anorm_drhoa)
2037 8586461 : e_ndr(ii) = e_ndr(ii) + 0.5_dp*exc_norm_drhoa
2038 : END IF
2039 :
2040 8586461 : s_bnorm_drhob = t25
2041 : u_x_bnorm_drhob = 0.2e1_dp*t370*t29*s_bnorm_drhob - 0.2e1_dp &
2042 8586461 : *t374*t376*s_bnorm_drhob
2043 8586461 : gx_bnorm_drhob = u_x_bnorm_drhob*t31 + u_x_b*u_x_bnorm_drhob*c_x_2
2044 8586461 : s_b_2norm_drhob = 0.2e1_dp*s_b*s_bnorm_drhob
2045 8586461 : s_avg_2norm_drhob = s_b_2norm_drhob/0.2e1_dp
2046 8586461 : t539 = t339*s_avg_2norm_drhob
2047 8586461 : u_c_abnorm_drhob = gamma_c_ab*s_avg_2norm_drhob*t162 - t337*t539
2048 8586461 : t543 = t486*s_b_2norm_drhob
2049 8586461 : 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 8586461 : u_c_abnorm_drhob*c_cab_2
2052 : gc_bnorm_drhob = u_c_bnorm_drhob*t176 + u_c_b*u_c_bnorm_drhob &
2053 8586461 : *c_css_2
2054 :
2055 8586461 : 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 8586461 : (e_lsda_c_ab*gc_abnorm_drhob + e_lsda_c_b*gc_bnorm_drhob)
2058 8586461 : e_ndr(ii) = e_ndr(ii) + 0.5_dp*exc_norm_drhob
2059 : END IF
2060 :
2061 8586461 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
2062 0 : t555 = t7**2
2063 0 : t560 = t186*my_rhoa
2064 0 : s_arhoarhoa = 0.28e2_dp/0.9e1_dp*my_norm_drhoa/t7/t560
2065 0 : t564 = s_arhoa**2
2066 0 : t568 = t194*s_a_2
2067 0 : t575 = t194*gamma_x
2068 0 : t576 = s_a_2**2
2069 0 : t577 = t575*t576
2070 0 : 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 0 : t577*t579*t564 - 0.2e1_dp*t196*t198*s_arhoarhoa
2074 0 : u_x_a1rhoa = u_x_arhoa
2075 0 : t600 = 0.1e1_dp/t207/rho
2076 0 : t601 = t33*t600
2077 0 : chirhoarhoa = -0.2e1_dp*t208 + 0.2e1_dp*t601
2078 0 : t605 = 0.3141592654e1_dp**2
2079 0 : t606 = 0.1e1_dp/t605
2080 0 : t608 = t207**2
2081 : rsrhoarhoa = -t4/t210/t36*t606/t608/0.18e2_dp + &
2082 0 : t4*t212*t600/0.6e1_dp
2083 0 : t619 = alpha_1_1*rsrhoa
2084 0 : t621 = t221*t235*t236
2085 0 : t626 = t40/t220/t51
2086 0 : t627 = t235**2
2087 0 : t631 = 0.1e1_dp/t46
2088 0 : t632 = beta_1_1*t631
2089 0 : t633 = rsrhoa**2
2090 0 : t639 = beta_3_1*t223
2091 0 : t644 = t48**2
2092 0 : t646 = rs**2
2093 0 : t647 = 0.1e1_dp/t646
2094 0 : t659 = t220**2
2095 0 : t661 = t40/t659
2096 0 : t662 = t55**2
2097 0 : 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 : 0.2e1_dp
2105 0 : e_c_u_01rhoa = e_c_u_0rhoa
2106 0 : t671 = alpha_1_2*rsrhoa
2107 0 : t673 = t244*t256*t257
2108 0 : t678 = t60/t243/t69
2109 0 : t679 = t256**2
2110 0 : t683 = beta_1_2*t631
2111 0 : t689 = beta_3_2*t223
2112 0 : t694 = t66**2
2113 0 : t707 = t243**2
2114 0 : t709 = t60/t707
2115 0 : t710 = t73**2
2116 0 : t711 = 0.1e1_dp/t710
2117 0 : t719 = alpha_1_3*rsrhoa
2118 0 : t721 = t265*t277*t278
2119 0 : t726 = t78/t264/t87
2120 0 : t727 = t277**2
2121 0 : t731 = beta_1_3*t631
2122 0 : t737 = beta_3_3*t223
2123 0 : t742 = t84**2
2124 0 : t755 = t264**2
2125 0 : t757 = t78/t755
2126 0 : t758 = t91**2
2127 0 : t759 = 0.1e1_dp/t758
2128 0 : alpha_c1rhoa = alpha_crhoa
2129 0 : t764 = t99**2
2130 0 : t765 = 0.1e1_dp/t764
2131 0 : t766 = chirhoa**2
2132 0 : t771 = t102**2
2133 0 : 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 : 0.4e1_dp/0.3e1_dp*t102*chirhoarhoa)*t97
2137 0 : f1rhoa = frhoa
2138 0 : t790 = alpha_c1rhoa*f
2139 0 : t793 = alpha_c*f1rhoa
2140 0 : t796 = t106*t107
2141 0 : t811 = e_c_u_1rhoa - e_c_u_01rhoa
2142 0 : t818 = t811*f
2143 0 : 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 0 : t766 + 0.4e1_dp*t113*t289*chirhoarhoa
2155 : epsilon_c_unif1rhoa = e_c_u_01rhoa + t790*t110 + t793*t110 - &
2156 0 : t293 + t818*t108 + t821*t108 + t301
2157 0 : t838 = t186**2
2158 : rs_arhoarhoa = -t4/t302/t116*t606/t838/0.18e2_dp + &
2159 0 : t4*t304/t560/0.6e1_dp
2160 0 : t858 = t327**2
2161 0 : t864 = rs_arhoa**2
2162 0 : t876 = rs_a**2
2163 0 : t877 = 0.1e1_dp/t876
2164 0 : t889 = t312**2
2165 0 : t892 = t133**2
2166 0 : epsilon_c_unif_a1rhoa = epsilon_c_unif_arhoa
2167 0 : s_a_2rhoarhoa = 0.2e1_dp*t564 + 0.2e1_dp*s_a*s_arhoarhoa
2168 0 : s_a_21rhoa = s_a_2rhoa
2169 0 : s_avg_2rhoarhoa = s_a_2rhoarhoa/0.2e1_dp
2170 0 : 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 0 : + epsilon_c_unif_a1rhoa
2180 0 : e_lsda_c_a1rhoa = epsilon_c_unif_a1rhoa*my_rhoa + epsilon_c_unif_a
2181 0 : t906 = t336*s_avg_2rhoa
2182 0 : t907 = t339*s_avg_21rhoa
2183 0 : t911 = t336*gamma_c_ab*s_avg_2
2184 0 : t913 = 0.1e1_dp/t338/t161
2185 0 : 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 0 : s_avg_2rhoarhoa
2189 0 : u_c_ab1rhoa = gamma_c_ab*s_avg_21rhoa*t162 - t337*t907
2190 0 : t925 = t344*s_a_2rhoa
2191 0 : t926 = t347*s_a_21rhoa
2192 0 : t929 = t344*gamma_c_ss
2193 0 : t930 = t929*s_a_2
2194 0 : t932 = 0.1e1_dp/t346/t164
2195 0 : 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 0 : s_a_2rhoarhoa
2199 0 : u_c_a1rhoa = gamma_c_ss*s_a_21rhoa*t165 - t345*t926
2200 0 : 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 0 : u_c_a1rhoa*c_css_2 + u_c_a*u_c_arhoarhoa*c_css_2))
2225 0 : e_r_r(ii) = e_r_r(ii) + 0.5_dp*0.5_dp*exc_rhoa_rhoa
2226 : END IF
2227 :
2228 0 : chirhoarhob = 0.2e1_dp*t601
2229 0 : rsrhoarhob = rsrhoarhoa
2230 0 : t974 = t221*t396*t236
2231 0 : t976 = alpha_1_1*rsrhob
2232 0 : t981 = rsrhoa*rsrhob
2233 0 : 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 : 0.2e1_dp
2241 0 : t1012 = t244*t410*t257
2242 0 : t1014 = alpha_1_2*rsrhob
2243 0 : t1047 = t265*t424*t278
2244 0 : 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 0 : t97
2249 0 : 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 : 0.4e1_dp*t113*t289*chirhoarhob
2261 : u_c_abrhoarhob = -0.2e1_dp*t906*t480 + 0.2e1_dp*t911*t914 &
2262 0 : *s_avg_2rhob
2263 :
2264 0 : 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 0 : u_c_abrhob*c_cab_2 + u_c_ab*u_c_abrhoarhob*c_cab_2))
2279 0 : e_r_r(ii) = e_r_r(ii) + 0.5_dp*exc_rhoa_rhob
2280 : END IF
2281 :
2282 0 : t1152 = t20**2
2283 0 : t1157 = t365*my_rhob
2284 0 : s_brhobrhob = 0.28e2_dp/0.9e1_dp*my_norm_drhob/t20/t1157
2285 0 : t1161 = s_brhob**2
2286 0 : t1165 = t194*s_b_2
2287 0 : t1172 = s_b_2**2
2288 0 : t1173 = t575*t1172
2289 0 : 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 0 : s_brhobrhob
2294 0 : u_x_b1rhob = u_x_brhob
2295 0 : chirhobrhob = 0.2e1_dp*t208 + 0.2e1_dp*t601
2296 0 : rsrhobrhob = rsrhoarhob
2297 0 : t1201 = t396**2
2298 0 : 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 0 : t1201*t663*t42/0.2e1_dp
2306 0 : e_c_u_01rhob = e_c_u_0rhob
2307 0 : t1236 = t410**2
2308 0 : t1270 = t424**2
2309 0 : alpha_c1rhob = alpha_crhob
2310 0 : 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 : 0.4e1_dp/0.3e1_dp*t102*chirhobrhob)*t97
2314 0 : f1rhob = frhob
2315 0 : t1321 = alpha_c1rhob*f
2316 0 : t1324 = alpha_c*f1rhob
2317 0 : t1341 = e_c_u_1rhob - e_c_u_01rhob
2318 0 : t1348 = t1341*f
2319 0 : 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 0 : *t107*t1299 + 0.4e1_dp*t113*t289*chirhobrhob
2331 : epsilon_c_unif1rhob = e_c_u_01rhob + t1321*t110 + t1324*t110 - &
2332 0 : t437 + t1348*t108 + t1351*t108 + t445
2333 0 : t1368 = t365**2
2334 : rs_brhobrhob = -t4/t446/t138*t606/t1368/0.18e2_dp &
2335 0 : + t4*t448/t1157/0.6e1_dp
2336 0 : t1388 = t471**2
2337 0 : t1394 = rs_brhob**2
2338 0 : t1406 = rs_b**2
2339 0 : t1407 = 0.1e1_dp/t1406
2340 0 : t1419 = t456**2
2341 0 : t1422 = t155**2
2342 0 : epsilon_c_unif_b1rhob = epsilon_c_unif_brhob
2343 0 : s_b_2rhobrhob = 0.2e1_dp*t1161 + 0.2e1_dp*s_b*s_brhobrhob
2344 0 : s_b_21rhob = s_b_2rhob
2345 0 : s_avg_2rhobrhob = s_b_2rhobrhob/0.2e1_dp
2346 0 : 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 0 : my_rhob + epsilon_c_unif_brhob + epsilon_c_unif_b1rhob
2356 0 : e_lsda_c_b1rhob = epsilon_c_unif_b1rhob*my_rhob + epsilon_c_unif_b
2357 0 : t1436 = t336*s_avg_2rhob
2358 0 : t1437 = t339*s_avg_21rhob
2359 0 : 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 0 : *s_avg_2rhobrhob
2363 0 : u_c_ab1rhob = gamma_c_ab*s_avg_21rhob*t162 - t337*t1437
2364 0 : t1451 = t344*s_b_2rhob
2365 0 : t1452 = t486*s_b_21rhob
2366 0 : t1455 = t929*s_b_2
2367 0 : t1457 = 0.1e1_dp/t485/t167
2368 0 : 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 0 : *s_b_2rhobrhob
2372 0 : u_c_b1rhob = gamma_c_ss*s_b_21rhob*t168 - t484*t1452
2373 :
2374 0 : 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 0 : *c_css_2))
2400 0 : e_r_r(ii) = e_r_r(ii) + 0.5_dp*0.5_dp*exc_rhob_rhob
2401 : END IF
2402 :
2403 0 : 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 : - 0.2e1_dp*t196*t198*s_arhoanorm_drhoa
2408 : s_a_2rhoanorm_drhoa = 0.2e1_dp*s_anorm_drhoa*s_arhoa + &
2409 0 : 0.2e1_dp*s_a*s_arhoanorm_drhoa
2410 0 : 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 0 : - 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 0 : t345*t347*s_a_2rhoanorm_drhoa
2417 :
2418 0 : 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 0 : *c_css_2 + u_c_a*u_c_arhoanorm_drhoa*c_css_2))
2428 0 : 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 0 : *t1440*s_avg_2norm_drhoa
2433 :
2434 0 : 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 0 : u_c_abrhobnorm_drhoa*c_cab_2))
2439 0 : e_r_ndr(ii) = e_r_ndr(ii) + 0.5_dp*0.5_dp*exc_rhob_norm_drhoa
2440 : END IF
2441 :
2442 0 : t1571 = s_anorm_drhoa**2
2443 : u_x_anorm_drhoanorm_drhoa = 0.2e1_dp*gamma_x*t1571*t16 - &
2444 0 : 0.10e2_dp*t568*t198*t1571 + 0.8e1_dp*t577*t579*t1571
2445 0 : s_a_2norm_drhoanorm_drhoa = 0.2e1_dp*t1571
2446 0 : s_a_21norm_drhoa = s_a_2norm_drhoa
2447 0 : s_avg_2norm_drhoanorm_drhoa = s_a_2norm_drhoanorm_drhoa/0.2e1_dp
2448 0 : s_avg_21norm_drhoa = s_a_21norm_drhoa/0.2e1_dp
2449 0 : t1589 = t336*s_avg_2norm_drhoa
2450 0 : t1590 = t339*s_avg_21norm_drhoa
2451 0 : 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 0 : s_avg_2norm_drhoanorm_drhoa
2456 0 : 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 0 : s_a_2norm_drhoanorm_drhoa
2461 :
2462 0 : 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 0 : t1605)*c_css_2 + u_c_a*u_c_anorm_drhoanorm_drhoa*c_css_2))
2472 0 : 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 0 : t914*s_avg_2norm_drhob
2477 :
2478 0 : 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 0 : u_c_abrhoanorm_drhob*c_cab_2))
2483 0 : e_r_ndr(ii) = e_r_ndr(ii) + 0.5_dp*0.5_dp*exc_rhoa_norm_drhob
2484 : END IF
2485 :
2486 0 : 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 0 : 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 : 0.2e1_dp*s_b*s_brhobnorm_drhob
2493 0 : 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 0 : 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 0 : - t484*t486*s_b_2rhobnorm_drhob
2500 :
2501 0 : 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 0 : *c_css_2 + u_c_b*u_c_brhobnorm_drhob*c_css_2))
2511 0 : 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 0 : t911*t1593*s_avg_2norm_drhob
2516 :
2517 0 : 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 0 : c_cab_2)
2522 0 : e_ndr_ndr(ii) = e_ndr_ndr(ii) + 0.5_dp*exc_norm_drhoa_norm_drhob
2523 : END IF
2524 :
2525 0 : t1719 = s_bnorm_drhob**2
2526 : u_x_bnorm_drhobnorm_drhob = 0.2e1_dp*gamma_x*t1719*t29 - &
2527 0 : 0.10e2_dp*t1165*t376*t1719 + 0.8e1_dp*t1173*t1175*t1719
2528 0 : s_b_2norm_drhobnorm_drhob = 0.2e1_dp*t1719
2529 0 : s_b_21norm_drhob = s_b_2norm_drhob
2530 0 : s_avg_2norm_drhobnorm_drhob = s_b_2norm_drhobnorm_drhob/0.2e1_dp
2531 0 : s_avg_21norm_drhob = s_b_21norm_drhob/0.2e1_dp
2532 0 : 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 0 : s_avg_2norm_drhobnorm_drhob
2538 0 : 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 0 : s_b_2norm_drhobnorm_drhob
2543 :
2544 0 : 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 0 : t1753)*c_css_2 + u_c_b*u_c_bnorm_drhobnorm_drhob*c_css_2))
2554 0 : 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 607 : END SUBROUTINE b97_lda_calc
2563 :
2564 : END MODULE xc_b97
|
