Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculates the exchange energy based on the Becke-Roussel exchange
10 : !> hole. Takes advantage of an analytical representation of the hole
11 : !> in order to avoid solving a non-linear equation by means of Newton-
12 : !> Raphson algorithm
13 : !> \note
14 : !> \par History
15 : !> 12.2008 created [mguidon]
16 : !> \author mguidon
17 : ! **************************************************************************************************
18 : MODULE xc_xbecke_roussel
19 : USE bibliography, ONLY: BeckeRoussel1989,&
20 : Proynov2007,&
21 : cite_reference
22 : USE input_section_types, ONLY: section_vals_type,&
23 : section_vals_val_get
24 : USE kinds, ONLY: dp
25 : USE mathconstants, ONLY: pi
26 : USE xc_derivative_desc, ONLY: &
27 : deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, deriv_norm_drho, &
28 : deriv_norm_drhoa, deriv_norm_drhob, deriv_rho, deriv_rhoa, deriv_rhob, deriv_tau, &
29 : deriv_tau_a, deriv_tau_b
30 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
31 : xc_dset_get_derivative
32 : USE xc_derivative_types, ONLY: xc_derivative_get,&
33 : xc_derivative_type
34 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
35 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
36 : xc_rho_set_type
37 : #include "../base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 : PRIVATE
41 :
42 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbecke_roussel'
44 :
45 : REAL(dp), PARAMETER, PRIVATE :: br_a1 = 1.5255251812009530_dp, &
46 : br_a2 = 0.4576575543602858_dp, &
47 : br_a3 = 0.4292036732051034_dp, &
48 : br_c0 = 0.7566445420735584_dp, &
49 : br_c1 = -2.6363977871370960_dp, &
50 : br_c2 = 5.4745159964232880_dp, &
51 : br_c3 = -12.657308127108290_dp, &
52 : br_c4 = 4.1250584725121360_dp, &
53 : br_c5 = -30.425133957163840_dp, &
54 : br_b0 = 0.4771976183772063_dp, &
55 : br_b1 = -1.7799813494556270_dp, &
56 : br_b2 = 3.8433841862302150_dp, &
57 : br_b3 = -9.5912050880518490_dp, &
58 : br_b4 = 2.1730180285916720_dp, &
59 : br_b5 = -30.425133851603660_dp, &
60 : br_d0 = 0.00004435009886795587_dp, &
61 : br_d1 = 0.58128653604457910_dp, &
62 : br_d2 = 66.742764515940610_dp, &
63 : br_d3 = 434.26780897229770_dp, &
64 : br_d4 = 824.7765766052239000_dp, &
65 : br_d5 = 1657.9652731582120_dp, &
66 : br_e0 = 0.00003347285060926091_dp, &
67 : br_e1 = 0.47917931023971350_dp, &
68 : br_e2 = 62.392268338574240_dp, &
69 : br_e3 = 463.14816427938120_dp, &
70 : br_e4 = 785.2360350104029000_dp, &
71 : br_e5 = 1657.962968223273000000_dp, &
72 : br_BB = 2.085749716493756_dp
73 :
74 : PUBLIC :: xbecke_roussel_lda_info, xbecke_roussel_lsd_info, xbecke_roussel_lda_eval, &
75 : xbecke_roussel_lsd_eval, x_br_lsd_y_lte_0, x_br_lsd_y_gt_0, x_br_lsd_y_lte_0_cutoff, &
76 : x_br_lsd_y_gt_0_cutoff
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief return various information on the functional
81 : !> \param reference string with the reference of the actual functional
82 : !> \param shortform string with the shortform of the functional name
83 : !> \param needs the components needed by this functional are set to
84 : !> true (does not set the unneeded components to false)
85 : !> \param max_deriv controls the number of derivatives
86 : !> \par History
87 : !> 12.2008 created [mguidon]
88 : !> \author mguidon
89 : ! **************************************************************************************************
90 99 : SUBROUTINE xbecke_roussel_lda_info(reference, shortform, needs, max_deriv)
91 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
92 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
93 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
94 :
95 99 : CALL cite_reference(BeckeRoussel1989)
96 99 : CALL cite_reference(Proynov2007)
97 99 : IF (PRESENT(reference)) THEN
98 : reference = "A.D. Becke, M.R. Roussel, "// &
99 3 : "Phys. Rev. A, vol. 39, n 8, pp. 3761-3767, (1989) {spin unpolarized}"
100 : END IF
101 99 : IF (PRESENT(shortform)) THEN
102 3 : shortform = "Becke-Roussel exchange hole (spin unpolarized)"
103 : END IF
104 99 : IF (PRESENT(needs)) THEN
105 96 : needs%rho = .TRUE.
106 96 : needs%norm_drho = .TRUE.
107 96 : needs%tau = .TRUE.
108 96 : needs%laplace_rho = .TRUE.
109 : END IF
110 :
111 99 : IF (PRESENT(max_deriv)) max_deriv = 1
112 :
113 99 : END SUBROUTINE xbecke_roussel_lda_info
114 :
115 : ! **************************************************************************************************
116 : !> \brief return various information on the functional
117 : !> \param reference string with the reference of the actual functional
118 : !> \param shortform string with the shortform of the functional name
119 : !> \param needs the components needed by this functional are set to
120 : !> true (does not set the unneeded components to false)
121 : !> \param max_deriv ...
122 : !> \par History
123 : !> 12.2008 created [mguidon]
124 : !> \author mguidon
125 : ! **************************************************************************************************
126 92 : SUBROUTINE xbecke_roussel_lsd_info(reference, shortform, needs, max_deriv)
127 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
128 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
129 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
130 :
131 92 : CALL cite_reference(BeckeRoussel1989)
132 92 : CALL cite_reference(Proynov2007)
133 92 : IF (PRESENT(reference)) THEN
134 : reference = "A.D. Becke, M.R. Roussel, "// &
135 2 : "Phys. Rev. A, vol. 39, n 8, pp. 3761-3767, (1989) {spin polarized}"
136 : END IF
137 92 : IF (PRESENT(shortform)) THEN
138 2 : shortform = "Becke-Roussel exchange hole (spin polarized)"
139 : END IF
140 92 : IF (PRESENT(needs)) THEN
141 90 : needs%rho_spin = .TRUE.
142 90 : needs%norm_drho_spin = .TRUE.
143 90 : needs%tau_spin = .TRUE.
144 90 : needs%laplace_rho_spin = .TRUE.
145 : END IF
146 92 : IF (PRESENT(max_deriv)) max_deriv = 1
147 :
148 92 : END SUBROUTINE xbecke_roussel_lsd_info
149 :
150 : ! **************************************************************************************************
151 : !> \brief evaluates the Becke Roussel exchange functional for lda
152 : !> \param rho_set the density where you want to evaluate the functional
153 : !> \param deriv_set place where to store the functional derivatives (they are
154 : !> added to the derivatives)
155 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
156 : !> \param br_params parameters for the becke roussel functional
157 : !> \par History
158 : !> 12.2008 created [mguidon]
159 : !> \author mguidon
160 : ! **************************************************************************************************
161 420 : SUBROUTINE xbecke_roussel_lda_eval(rho_set, deriv_set, grad_deriv, br_params)
162 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
163 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
164 : INTEGER, INTENT(in) :: grad_deriv
165 : TYPE(section_vals_type), POINTER :: br_params
166 :
167 : CHARACTER(len=*), PARAMETER :: routineN = 'xbecke_roussel_lda_eval'
168 :
169 : INTEGER :: handle, npoints
170 : INTEGER, DIMENSION(2, 3) :: bo
171 : REAL(dp) :: gamma, R, sx
172 : REAL(kind=dp) :: epsilon_rho
173 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
174 84 : POINTER :: dummy, e_0, e_laplace_rho, e_ndrho, &
175 84 : e_rho, e_tau, laplace_rho, norm_drho, &
176 84 : rho, tau
177 : TYPE(xc_derivative_type), POINTER :: deriv
178 :
179 84 : CALL timeset(routineN, handle)
180 :
181 : CALL xc_rho_set_get(rho_set, rho=rho, norm_drho=norm_drho, &
182 : tau=tau, laplace_rho=laplace_rho, local_bounds=bo, &
183 84 : rho_cutoff=epsilon_rho)
184 84 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
185 :
186 84 : dummy => rho
187 :
188 84 : e_0 => dummy
189 84 : e_rho => dummy
190 84 : e_ndrho => dummy
191 84 : e_tau => dummy
192 84 : e_laplace_rho => dummy
193 :
194 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
195 84 : allocate_deriv=.TRUE.)
196 84 : CALL xc_derivative_get(deriv, deriv_data=e_0)
197 84 : IF (grad_deriv >= 1) THEN
198 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
199 84 : allocate_deriv=.TRUE.)
200 84 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
201 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
202 84 : allocate_deriv=.TRUE.)
203 84 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
204 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
205 84 : allocate_deriv=.TRUE.)
206 84 : CALL xc_derivative_get(deriv, deriv_data=e_tau)
207 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho], &
208 84 : allocate_deriv=.TRUE.)
209 84 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho)
210 : END IF
211 84 : IF (grad_deriv > 1) THEN
212 0 : CPABORT("derivatives bigger than 1 not implemented")
213 : END IF
214 :
215 84 : CALL section_vals_val_get(br_params, "scale_x", r_val=sx)
216 84 : CALL section_vals_val_get(br_params, "CUTOFF_RADIUS", r_val=R)
217 84 : CALL section_vals_val_get(br_params, "GAMMA", r_val=gamma)
218 :
219 : !$OMP PARALLEL DEFAULT(NONE) &
220 : !$OMP SHARED(rho, norm_drho, laplace_rho, tau, e_0, e_rho) &
221 : !$OMP SHARED(e_ndrho, e_tau, e_laplace_rho, grad_deriv) &
222 : !$OMP SHARED(npoints, epsilon_rho) &
223 84 : !$OMP SHARED(sx, r, gamma)
224 :
225 : CALL xbecke_roussel_lda_calc(rho=rho, norm_drho=norm_drho, &
226 : laplace_rho=laplace_rho, tau=tau, e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
227 : e_tau=e_tau, e_laplace_rho=e_laplace_rho, grad_deriv=grad_deriv, &
228 : npoints=npoints, epsilon_rho=epsilon_rho, &
229 : sx=sx, R=R, gamma=gamma)
230 :
231 : !$OMP END PARALLEL
232 :
233 84 : CALL timestop(handle)
234 84 : END SUBROUTINE xbecke_roussel_lda_eval
235 :
236 : ! **************************************************************************************************
237 : !> \brief Precalculates which branch of the code has to be taken
238 : !> There are two main branches, one for a truncated potential and a cutoff
239 : !> radius, the other for the full coulomb interaction. In the end, the code
240 : !> for the lsd part will be called and the lda part is obtained via spin
241 : !> scaling relations
242 : !> \param rho grid values
243 : !> \param norm_drho grid values
244 : !> \param laplace_rho grid values
245 : !> \param tau grid values
246 : !> \param e_0 energies and derivatives
247 : !> \param e_rho energies and derivatives
248 : !> \param e_ndrho energies and derivatives
249 : !> \param e_tau energies and derivatives
250 : !> \param e_laplace_rho energies and derivatives
251 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
252 : !> \param npoints size of the grids
253 : !> \param epsilon_rho cutoffs
254 : !> \param sx scales the exchange potential and energies
255 : !> \param R cutoff Radius for truncated case
256 : !> \param gamma parameter from original publication, usually set to 1
257 : !> \par History
258 : !> 12.2008 created [mguidon]
259 : !> \author mguidon
260 : ! **************************************************************************************************
261 84 : SUBROUTINE xbecke_roussel_lda_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
262 84 : e_ndrho, e_tau, e_laplace_rho, grad_deriv, npoints, &
263 : epsilon_rho, sx, R, gamma)
264 :
265 : INTEGER, INTENT(in) :: npoints, grad_deriv
266 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_laplace_rho, e_tau, e_ndrho, e_rho, e_0
267 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: tau, laplace_rho, norm_drho, rho
268 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R, gamma
269 :
270 : INTEGER :: ip
271 : REAL(dp) :: e_diff, e_old, my_laplace_rho, my_ndrho, &
272 : my_rho, my_tau, t1, t15, t16, t2, t3, &
273 : t4, t5, t8, t9, yval
274 :
275 : ! Precalculate y, in order to chose the correct branch afterwards
276 :
277 : !$OMP DO
278 :
279 : DO ip = 1, npoints
280 1283968 : my_rho = 0.5_dp*MAX(rho(ip), 0.0_dp)
281 1283968 : IF (my_rho > epsilon_rho) THEN
282 1283968 : my_ndrho = 0.5_dp*MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp)
283 1283968 : my_tau = 0.5_dp*MAX(EPSILON(0.0_dp)*1.e4_dp, tau(ip))
284 1283968 : my_laplace_rho = 0.5_dp*laplace_rho(ip)
285 :
286 1283968 : t1 = pi**(0.1e1_dp/0.3e1_dp)
287 1283968 : t2 = t1**2
288 1283968 : t3 = my_rho**(0.1e1_dp/0.3e1_dp)
289 1283968 : t4 = t3**2
290 1283968 : t5 = t4*my_rho
291 1283968 : t8 = my_ndrho**2
292 1283968 : t9 = 0.1e1_dp/my_rho
293 : ! *** CP2K defines tau in a different way as compared to Becke !!!
294 1283968 : t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
295 1283968 : t16 = 0.1e1_dp/t15
296 1283968 : yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
297 1283968 : IF (R == 0.0_dp) THEN
298 890752 : IF (yval <= 0.0_dp) THEN
299 138801 : e_old = e_0(ip)
300 : CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), &
301 : e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), &
302 138801 : sx, gamma, grad_deriv)
303 : ! VERY UGLY HACK e_0 has to multiplied by the factor 2
304 138801 : e_diff = e_0(ip) - e_old
305 138801 : e_0(ip) = e_0(ip) + e_diff
306 : ELSE
307 751951 : e_old = e_0(ip)
308 : CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), &
309 : e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), &
310 751951 : sx, gamma, grad_deriv)
311 : ! VERY UGLY HACK e_0 has to multiplied by the factor 2
312 751951 : e_diff = e_0(ip) - e_old
313 751951 : e_0(ip) = e_0(ip) + e_diff
314 : END IF
315 : ELSE
316 393216 : IF (yval <= 0.0_dp) THEN
317 10559 : e_old = e_0(ip)
318 : CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), &
319 : e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), &
320 10559 : sx, R, gamma, grad_deriv)
321 : ! VERY UGLY HACK e_0 has to multiplied by the factor 2
322 10559 : e_diff = e_0(ip) - e_old
323 10559 : e_0(ip) = e_0(ip) + e_diff
324 : ELSE
325 382657 : e_old = e_0(ip)
326 : CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), &
327 : e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), &
328 382657 : sx, R, gamma, grad_deriv)
329 : ! VERY UGLY HACK e_0 has to multiplied by the factor 2
330 382657 : e_diff = e_0(ip) - e_old
331 382657 : e_0(ip) = e_0(ip) + e_diff
332 : END IF
333 : END IF
334 : END IF
335 : END DO
336 :
337 : !$OMP END DO
338 :
339 84 : END SUBROUTINE xbecke_roussel_lda_calc
340 :
341 : ! **************************************************************************************************
342 : !> \brief evaluates the Becke Roussel exchange functional for lda
343 : !> \param rho_set the density where you want to evaluate the functional
344 : !> \param deriv_set place where to store the functional derivatives (they are
345 : !> added to the derivatives)
346 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
347 : !> \param br_params parameters for the becke roussel functional
348 : !> \par History
349 : !> 12.2008 created [mguidon]
350 : !> \author mguidon
351 : ! **************************************************************************************************
352 410 : SUBROUTINE xbecke_roussel_lsd_eval(rho_set, deriv_set, grad_deriv, br_params)
353 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
354 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
355 : INTEGER, INTENT(in) :: grad_deriv
356 : TYPE(section_vals_type), POINTER :: br_params
357 :
358 : CHARACTER(len=*), PARAMETER :: routineN = 'xbecke_roussel_lsd_eval'
359 :
360 : INTEGER :: handle, npoints
361 : INTEGER, DIMENSION(2, 3) :: bo
362 : REAL(dp) :: gamma, R, sx
363 : REAL(kind=dp) :: epsilon_rho
364 82 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, &
365 82 : e_laplace_rhob, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, e_tau_a, e_tau_b, laplace_rhoa, &
366 82 : laplace_rhob, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
367 : TYPE(xc_derivative_type), POINTER :: deriv
368 :
369 82 : CALL timeset(routineN, handle)
370 :
371 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
372 : norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, laplace_rhoa=laplace_rhoa, &
373 : laplace_rhob=laplace_rhob, local_bounds=bo, &
374 82 : rho_cutoff=epsilon_rho)
375 82 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
376 :
377 82 : dummy => rhoa
378 :
379 82 : e_0 => dummy
380 82 : e_rhoa => dummy
381 82 : e_rhob => dummy
382 82 : e_ndrhoa => dummy
383 82 : e_ndrhob => dummy
384 82 : e_tau_a => dummy
385 82 : e_tau_b => dummy
386 82 : e_laplace_rhoa => dummy
387 82 : e_laplace_rhob => dummy
388 :
389 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
390 82 : allocate_deriv=.TRUE.)
391 82 : CALL xc_derivative_get(deriv, deriv_data=e_0)
392 82 : IF (grad_deriv >= 1) THEN
393 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
394 54 : allocate_deriv=.TRUE.)
395 54 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
396 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
397 54 : allocate_deriv=.TRUE.)
398 54 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
399 :
400 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
401 54 : allocate_deriv=.TRUE.)
402 54 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
403 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
404 54 : allocate_deriv=.TRUE.)
405 54 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
406 :
407 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], &
408 54 : allocate_deriv=.TRUE.)
409 54 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a)
410 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], &
411 54 : allocate_deriv=.TRUE.)
412 54 : CALL xc_derivative_get(deriv, deriv_data=e_tau_b)
413 :
414 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa], &
415 54 : allocate_deriv=.TRUE.)
416 54 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa)
417 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob], &
418 54 : allocate_deriv=.TRUE.)
419 54 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob)
420 : END IF
421 82 : IF (grad_deriv > 1) THEN
422 0 : CPABORT("derivatives bigger than 1 not implemented")
423 : END IF
424 :
425 82 : CALL section_vals_val_get(br_params, "scale_x", r_val=sx)
426 82 : CALL section_vals_val_get(br_params, "CUTOFF_RADIUS", r_val=R)
427 82 : CALL section_vals_val_get(br_params, "GAMMA", r_val=gamma)
428 :
429 : !$OMP PARALLEL DEFAULT (NONE) &
430 : !$OMP SHARED(rhoa, norm_drhoa, laplace_rhoa, tau_a, e_0) &
431 : !$OMP SHARED(e_rhoa, e_ndrhoa, e_tau_a, e_laplace_rhoa) &
432 : !$OMP SHARED(grad_deriv, npoints, epsilon_rho) &
433 : !$OMP SHARED(sx, r, gamma) &
434 : !$OMP SHARED(rhob, norm_drhob, laplace_rhob, tau_b, e_rhob) &
435 82 : !$OMP SHARED(e_ndrhob, e_tau_b, e_laplace_rhob)
436 :
437 : CALL xbecke_roussel_lsd_calc(rho=rhoa, norm_drho=norm_drhoa, &
438 : laplace_rho=laplace_rhoa, tau=tau_a, e_0=e_0, e_rho=e_rhoa, e_ndrho=e_ndrhoa, &
439 : e_tau=e_tau_a, e_laplace_rho=e_laplace_rhoa, grad_deriv=grad_deriv, &
440 : npoints=npoints, epsilon_rho=epsilon_rho, &
441 : sx=sx, R=R, gamma=gamma)
442 :
443 : CALL xbecke_roussel_lsd_calc(rho=rhob, norm_drho=norm_drhob, &
444 : laplace_rho=laplace_rhob, tau=tau_b, e_0=e_0, e_rho=e_rhob, e_ndrho=e_ndrhob, &
445 : e_tau=e_tau_b, e_laplace_rho=e_laplace_rhob, grad_deriv=grad_deriv, &
446 : npoints=npoints, epsilon_rho=epsilon_rho, &
447 : sx=sx, R=R, gamma=gamma)
448 :
449 : !$OMP END PARALLEL
450 :
451 82 : CALL timestop(handle)
452 82 : END SUBROUTINE xbecke_roussel_lsd_eval
453 :
454 : ! **************************************************************************************************
455 : !> \brief Precalculates which branch of the code has to be taken
456 : !> There are two main branches, one for a truncated potential and a cutoff
457 : !> radius, the other for the full coulomb interaction
458 : !> \param rho grid values
459 : !> \param norm_drho grid values
460 : !> \param laplace_rho grid values
461 : !> \param tau grid values
462 : !> \param e_0 energies and derivatives
463 : !> \param e_rho energies and derivatives
464 : !> \param e_ndrho energies and derivatives
465 : !> \param e_tau energies and derivatives
466 : !> \param e_laplace_rho energies and derivatives
467 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
468 : !> \param npoints size of the grids
469 : !> \param epsilon_rho cutoffs
470 : !> \param sx scales the exchange potential and energies
471 : !> \param R cutoff Radius for truncated case
472 : !> \param gamma parameter from original publication, usually set to 1
473 : !> \par History
474 : !> 12.2008 created [mguidon]
475 : !> \author mguidon
476 : ! **************************************************************************************************
477 164 : SUBROUTINE xbecke_roussel_lsd_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
478 164 : e_ndrho, e_tau, e_laplace_rho, grad_deriv, npoints, &
479 : epsilon_rho, sx, R, gamma)
480 :
481 : INTEGER, INTENT(in) :: npoints, grad_deriv
482 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_laplace_rho, e_tau, e_ndrho, e_rho, e_0
483 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: tau, laplace_rho, norm_drho, rho
484 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R, gamma
485 :
486 : INTEGER :: ip
487 : REAL(dp) :: my_laplace_rho, my_ndrho, my_rho, &
488 : my_tau, t1, t15, t16, t2, t3, t4, t5, &
489 : t8, t9, yval
490 :
491 : ! Precalculate y, in order to chose the correct branch afterwards
492 :
493 : !$OMP DO
494 :
495 : DO ip = 1, npoints
496 7472250 : my_rho = MAX(rho(ip), 0.0_dp)
497 7472250 : IF (my_rho > epsilon_rho) THEN
498 7472244 : my_ndrho = MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp)
499 7472244 : my_tau = 1.0_dp*MAX(EPSILON(0.0_dp)*1.e4_dp, tau(ip))
500 7472244 : my_laplace_rho = 1.0_dp*laplace_rho(ip)
501 :
502 7472244 : t1 = pi**(0.1e1_dp/0.3e1_dp)
503 7472244 : t2 = t1**2
504 7472244 : t3 = my_rho**(0.1e1_dp/0.3e1_dp)
505 7472244 : t4 = t3**2
506 7472244 : t5 = t4*my_rho
507 7472244 : t8 = my_ndrho**2
508 7472244 : t9 = 0.1e1_dp/my_rho
509 : ! *** CP2K defines tau in a different way as compared to Becke !!!
510 7472244 : t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
511 7472244 : t16 = 0.1e1_dp/t15
512 7472244 : yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
513 7472244 : IF (R == 0.0_dp) THEN
514 2187000 : IF (yval <= 0.0_dp) THEN
515 : CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), &
516 : e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), &
517 48750 : sx, gamma, grad_deriv)
518 : ELSE
519 : CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), &
520 : e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), &
521 2138250 : sx, gamma, grad_deriv)
522 : END IF
523 : ELSE
524 5285244 : IF (yval <= 0.0_dp) THEN
525 : CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), &
526 : e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), &
527 117314 : sx, R, gamma, grad_deriv)
528 : ELSE
529 : CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), &
530 : e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), &
531 5167930 : sx, R, gamma, grad_deriv)
532 : END IF
533 : END IF
534 : END IF
535 : END DO
536 :
537 : !$OMP END DO
538 :
539 164 : END SUBROUTINE xbecke_roussel_lsd_calc
540 :
541 : ! **************************************************************************************************
542 : !> \brief full range evaluation for y <= 0
543 : !> \param rho ...
544 : !> \param ndrho ...
545 : !> \param tau ...
546 : !> \param laplace_rho ...
547 : !> \param e_0 ...
548 : !> \param e_rho ...
549 : !> \param e_ndrho ...
550 : !> \param e_tau ...
551 : !> \param e_laplace_rho ...
552 : !> \param sx ...
553 : !> \param gamma ...
554 : !> \param grad_deriv ...
555 : !> \par History
556 : !> 12.2008 created [mguidon]
557 : !> \author mguidon
558 : ! **************************************************************************************************
559 187551 : SUBROUTINE x_br_lsd_y_lte_0(rho, ndrho, tau, laplace_rho, e_0, &
560 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
561 : sx, gamma, grad_deriv)
562 : REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho
563 : REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho
564 : REAL(dp), INTENT(IN) :: sx, gamma
565 : INTEGER, INTENT(IN) :: grad_deriv
566 :
567 : REAL(KIND=dp) :: t1, t100, t101, t102, t108, t111, t113, t114, t117, t118, t120, t121, t122, &
568 : t129, t130, t134, t135, t138, t142, t143, t146, t147, t150, t152, t153, t157, t158, t16, &
569 : t161, t164, t165, t166, t169, t17, t170, t172, t173, t19, t199, t2, t20, t202, t207, &
570 : t208, t209, t21, t217, t218, t22, t220, t227, t229, t23, t234, t246, t249, t252, t255, &
571 : t259, t26, t263, t267, t27, t271, t274, t275, t276, t28, t29, t293, t295, t3, t304, t307, &
572 : t308, t32, t320, t33, t34, t340, t341, t342, t344, t346, t349, t35, t350, t353, t354, &
573 : t357, t358, t36, t361, t362, t365, t366, t367, t37, t379, t38
574 : REAL(KIND=dp) :: t381, t387, t39, t4, t401, t42, t422, t43, t434, t435, t436, t44, t448, &
575 : t45, t450, t47, t471, t48, t5, t51, t52, t53, t54, t55, t56, t57, t61, t62, t63, t64, &
576 : t66, t67, t70, t71, t72, t75, t78, t81, t84, t87, t88, t89, t9, t91, t92, t93, t94, t95, &
577 : t96, t97, yval
578 :
579 187551 : t1 = pi**(0.1e1_dp/0.3e1_dp)
580 187551 : t2 = t1**2
581 187551 : t3 = rho**(0.1e1_dp/0.3e1_dp)
582 187551 : t4 = t3**2
583 187551 : t5 = t4*rho
584 187551 : t9 = ndrho**2
585 187551 : t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9/rho/0.4e1_dp)/0.3e1_dp
586 187551 : t17 = 0.1e1_dp/t16
587 187551 : yval = 0.2e1_dp/0.3e1_dp*t2*t5*t17
588 187551 : t19 = t3*rho
589 187551 : t20 = 0.3141592654e1_dp**(0.1e1_dp/0.3e1_dp)
590 187551 : t21 = t19*t20
591 187551 : t22 = br_a1*t2
592 187551 : t23 = t5*t17
593 187551 : t26 = 0.2e1_dp/0.3e1_dp*t22*t23 + br_a2
594 187551 : t27 = ATAN(t26)
595 187551 : t28 = -t27 + br_a3
596 187551 : t29 = br_c1*t2
597 187551 : t32 = t1*pi
598 187551 : t33 = br_c2*t32
599 187551 : t34 = rho**2
600 187551 : t35 = t34*rho
601 187551 : t36 = t3*t35
602 187551 : t37 = t16**2
603 187551 : t38 = 0.1e1_dp/t37
604 187551 : t39 = t36*t38
605 187551 : t42 = pi**2
606 187551 : t43 = br_c3*t42
607 187551 : t44 = t34**2
608 187551 : t45 = t44*rho
609 187551 : t47 = 0.1e1_dp/t37/t16
610 187551 : t48 = t45*t47
611 187551 : t51 = t2*t42
612 187551 : t52 = br_c4*t51
613 187551 : t53 = t44*t34
614 187551 : t54 = t4*t53
615 187551 : t55 = t37**2
616 187551 : t56 = 0.1e1_dp/t55
617 187551 : t57 = t54*t56
618 187551 : t61 = t1*t42*pi
619 187551 : t62 = br_c5*t61
620 187551 : t63 = t44**2
621 187551 : t64 = t3*t63
622 187551 : t66 = 0.1e1_dp/t55/t16
623 187551 : t67 = t64*t66
624 : t70 = br_c0 + 0.2e1_dp/0.3e1_dp*t29*t23 + 0.4e1_dp/0.9e1_dp*t33*t39 &
625 : + 0.8e1_dp/0.27e2_dp*t43*t48 + 0.16e2_dp/0.81e2_dp*t52*t57 + 0.32e2_dp &
626 187551 : /0.243e3_dp*t62*t67
627 187551 : t71 = t28*t70
628 187551 : t72 = br_b1*t2
629 187551 : t75 = br_b2*t32
630 187551 : t78 = br_b3*t42
631 187551 : t81 = br_b4*t51
632 187551 : t84 = br_b5*t61
633 : t87 = br_b0 + 0.2e1_dp/0.3e1_dp*t72*t23 + 0.4e1_dp/0.9e1_dp*t75*t39 &
634 : + 0.8e1_dp/0.27e2_dp*t78*t48 + 0.16e2_dp/0.81e2_dp*t81*t57 + 0.32e2_dp &
635 187551 : /0.243e3_dp*t84*t67
636 187551 : t88 = 0.1e1_dp/t87
637 187551 : t89 = t71*t88
638 187551 : t91 = EXP(t89/0.3e1_dp)
639 187551 : t92 = t21*t91
640 187551 : t93 = 0.1e1_dp/t28
641 187551 : t94 = 0.1e1_dp/t70
642 187551 : t95 = t93*t94
643 187551 : t96 = EXP(-t89)
644 187551 : t97 = t88*t96
645 187551 : t100 = 0.1e1_dp - t96 - t71*t97/0.2e1_dp
646 187551 : t101 = t87*t100
647 187551 : t102 = t95*t101
648 187551 : e_0 = e_0 + (-t92*t102)*sx
649 187551 : IF (grad_deriv >= 1) THEN
650 187551 : t108 = t4*t17
651 187551 : t111 = 0.1e1_dp/t3
652 187551 : t113 = t38*gamma
653 187551 : t114 = t113*t9
654 187551 : t117 = 0.10e2_dp/0.9e1_dp*t22*t108 + t22*t111*t114/0.18e2_dp
655 187551 : t118 = t26**2
656 187551 : t120 = 0.1e1_dp/(0.1e1_dp + t118)
657 187551 : t121 = t117*t120
658 187551 : t122 = t70*t88
659 187551 : t129 = t3*t34
660 187551 : t130 = t129*t38
661 187551 : t134 = t47*gamma
662 187551 : t135 = t134*t9
663 187551 : t138 = t44*t47
664 187551 : t142 = t56*gamma
665 187551 : t143 = t142*t9
666 187551 : t146 = t4*t45
667 187551 : t147 = t146*t56
668 187551 : t150 = t4*t44
669 187551 : t152 = t66*gamma
670 187551 : t153 = t152*t9
671 187551 : t157 = t3*t44*t35
672 187551 : t158 = t157*t66
673 187551 : t161 = t3*t53
674 187551 : t164 = 0.1e1_dp/t55/t37
675 187551 : t165 = t164*gamma
676 187551 : t166 = t165*t9
677 : t169 = 0.10e2_dp/0.9e1_dp*t29*t108 + t29*t111*t114/0.18e2_dp + 0.40e2_dp &
678 : /0.27e2_dp*t33*t130 + 0.2e1_dp/0.27e2_dp*t33*t19*t135 + &
679 : 0.40e2_dp/0.27e2_dp*t43*t138 + 0.2e1_dp/0.27e2_dp*t43*t35*t143 + &
680 : 0.320e3_dp/0.243e3_dp*t52*t147 + 0.16e2_dp/0.243e3_dp*t52*t150* &
681 : t153 + 0.800e3_dp/0.729e3_dp*t62*t158 + 0.40e2_dp/0.729e3_dp*t62*t161 &
682 187551 : *t166
683 187551 : t170 = t28*t169
684 187551 : t172 = t87**2
685 187551 : t173 = 0.1e1_dp/t172
686 : t199 = 0.10e2_dp/0.9e1_dp*t72*t108 + t72*t111*t114/0.18e2_dp + 0.40e2_dp &
687 : /0.27e2_dp*t75*t130 + 0.2e1_dp/0.27e2_dp*t75*t19*t135 + &
688 : 0.40e2_dp/0.27e2_dp*t78*t138 + 0.2e1_dp/0.27e2_dp*t78*t35*t143 + &
689 : 0.320e3_dp/0.243e3_dp*t81*t147 + 0.16e2_dp/0.243e3_dp*t81*t150* &
690 : t153 + 0.800e3_dp/0.729e3_dp*t84*t158 + 0.40e2_dp/0.729e3_dp*t84*t161 &
691 187551 : *t166
692 187551 : t202 = -t121*t122 + t170*t88 - t71*t173*t199
693 187551 : t207 = t28**2
694 187551 : t208 = 0.1e1_dp/t207
695 187551 : t209 = t91*t208
696 187551 : t217 = t21*t91*t93
697 187551 : t218 = t70**2
698 187551 : t220 = 0.1e1_dp/t218*t87
699 187551 : t227 = -t202
700 187551 : t229 = t122*t96
701 187551 : t234 = t173*t96
702 : e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t3*t20*t91*t102 - t21*t202*t91* &
703 : t102/0.3e1_dp - t21*t209*t94*t87*t100*t117*t120 + t217 &
704 : *t220*t100*t169 - t92*t95*t199*t100 - t92*t95*t87* &
705 : (-t227*t96 + t121*t229/0.2e1_dp - t170*t97/0.2e1_dp + t71*t234 &
706 187551 : *t199/0.2e1_dp - t71*t88*t227*t96/0.2e1_dp))*sx
707 187551 : t246 = t4*t38
708 187551 : t249 = t120*t70
709 187551 : t252 = t22*t246*gamma*ndrho*t249*t88
710 187551 : t255 = t113*ndrho
711 187551 : t259 = t134*ndrho
712 187551 : t263 = t142*ndrho
713 187551 : t267 = t152*ndrho
714 187551 : t271 = t165*ndrho
715 : t274 = -t29*t4*t255/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t33*t129*t259 &
716 : - 0.4e1_dp/0.27e2_dp*t43*t44*t263 - 0.32e2_dp/0.243e3_dp*t52*t146 &
717 187551 : *t267 - 0.80e2_dp/0.729e3_dp*t62*t157*t271
718 187551 : t275 = t28*t274
719 187551 : t276 = t275*t88
720 : t293 = -t72*t4*t255/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t75*t129*t259 &
721 : - 0.4e1_dp/0.27e2_dp*t78*t44*t263 - 0.32e2_dp/0.243e3_dp*t81*t146 &
722 187551 : *t267 - 0.80e2_dp/0.729e3_dp*t84*t157*t271
723 187551 : t295 = t71*t173*t293
724 187551 : t304 = t208*t94*t87
725 187551 : t307 = t100*br_a1*t2
726 187551 : t308 = ndrho*t120
727 187551 : t320 = -t252/0.9e1_dp - t276 + t295
728 : e_ndrho = e_ndrho + (-t21*(t252/0.27e2_dp + t276/0.3e1_dp - t295/0.3e1_dp)*t91 &
729 : *t102 + t34*t20*t91*t304*t307*t113*t308/0.9e1_dp + t217 &
730 : *t220*t100*t274 - t92*t95*t293*t100 - t92*t95*t87 &
731 : *(-t320*t96 - t22*t246*gamma*t308*t229/0.18e2_dp - t275 &
732 : *t97/0.2e1_dp + t71*t234*t293/0.2e1_dp - t71*t88*t320*t96 &
733 187551 : /0.2e1_dp))*sx
734 187551 : t340 = t5*t38
735 187551 : t341 = t22*t340
736 187551 : t342 = gamma*t120
737 187551 : t344 = t341*t342*t122
738 187551 : t346 = t340*gamma
739 187551 : t349 = t36*t47
740 187551 : t350 = t349*gamma
741 187551 : t353 = t45*t56
742 187551 : t354 = t353*gamma
743 187551 : t357 = t54*t66
744 187551 : t358 = t357*gamma
745 187551 : t361 = t64*t164
746 187551 : t362 = t361*gamma
747 : t365 = 0.4e1_dp/0.9e1_dp*t29*t346 + 0.16e2_dp/0.27e2_dp*t33*t350 + &
748 : 0.16e2_dp/0.27e2_dp*t43*t354 + 0.128e3_dp/0.243e3_dp*t52*t358 + 0.320e3_dp &
749 187551 : /0.729e3_dp*t62*t362
750 187551 : t366 = t28*t365
751 187551 : t367 = t366*t88
752 : t379 = 0.4e1_dp/0.9e1_dp*t72*t346 + 0.16e2_dp/0.27e2_dp*t75*t350 + &
753 : 0.16e2_dp/0.27e2_dp*t78*t354 + 0.128e3_dp/0.243e3_dp*t81*t358 + 0.320e3_dp &
754 187551 : /0.729e3_dp*t84*t362
755 187551 : t381 = t71*t173*t379
756 187551 : t387 = t35*t20
757 187551 : t401 = 0.4e1_dp/0.9e1_dp*t344 - t367 + t381
758 : e_tau = e_tau + (-t21*(-0.4e1_dp/0.27e2_dp*t344 + t367/0.3e1_dp - t381/0.3e1_dp) &
759 : *t91*t102 - 0.4e1_dp/0.9e1_dp*t387*t91*t304*t307*t113* &
760 : t120 + t217*t220*t100*t365 - t92*t95*t379*t100 - t92 &
761 : *t95*t87*(-t401*t96 + 0.2e1_dp/0.9e1_dp*t341*t342*t229 - &
762 : t366*t97/0.2e1_dp + t71*t234*t379/0.2e1_dp - t71*t88*t401 &
763 187551 : *t96/0.2e1_dp))*sx
764 187551 : t422 = t22*t5*t38*t120*t122
765 : t434 = -t29*t340/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t33*t349 - 0.4e1_dp/ &
766 : 0.27e2_dp*t43*t353 - 0.32e2_dp/0.243e3_dp*t52*t357 - 0.80e2_dp/0.729e3_dp &
767 187551 : *t62*t361
768 187551 : t435 = t28*t434
769 187551 : t436 = t435*t88
770 : t448 = -t72*t340/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t75*t349 - 0.4e1_dp/ &
771 : 0.27e2_dp*t78*t353 - 0.32e2_dp/0.243e3_dp*t81*t357 - 0.80e2_dp/0.729e3_dp &
772 187551 : *t84*t361
773 187551 : t450 = t71*t173*t448
774 187551 : t471 = -t422/0.9e1_dp - t436 + t450
775 : e_laplace_rho = e_laplace_rho + (-t21*(t422/0.27e2_dp + t436/0.3e1_dp - t450/0.3e1_dp)*t91* &
776 : t102 + t387*t209*t94*t101*br_a1*t2*t38*t120/0.9e1_dp &
777 : + t217*t220*t100*t434 - t92*t95*t448*t100 - t92*t95 &
778 : *t87*(-t471*t96 - t341*t249*t97/0.18e2_dp - t435*t97/0.2e1_dp &
779 187551 : + t71*t234*t448/0.2e1_dp - t71*t88*t471*t96/0.2e1_dp))*sx
780 : END IF
781 187551 : END SUBROUTINE x_br_lsd_y_lte_0
782 :
783 : ! **************************************************************************************************
784 : !> \brief Full range evaluation for y > 0
785 : !> \param rho ...
786 : !> \param ndrho ...
787 : !> \param tau ...
788 : !> \param laplace_rho ...
789 : !> \param e_0 ...
790 : !> \param e_rho ...
791 : !> \param e_ndrho ...
792 : !> \param e_tau ...
793 : !> \param e_laplace_rho ...
794 : !> \param sx ...
795 : !> \param gamma ...
796 : !> \param grad_deriv ...
797 : !> \par History
798 : !> 12.2008 created [mguidon]
799 : !> \author mguidon
800 : ! **************************************************************************************************
801 2890201 : SUBROUTINE x_br_lsd_y_gt_0(rho, ndrho, tau, laplace_rho, e_0, &
802 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
803 : sx, gamma, grad_deriv)
804 : REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho
805 : REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho
806 : REAL(dp), INTENT(IN) :: sx, gamma
807 : INTEGER, INTENT(IN) :: grad_deriv
808 :
809 : REAL(KIND=dp) :: t1, t102, t103, t104, t106, t107, t108, t109, t110, t111, t112, t115, t117, &
810 : t124, t131, t134, t137, t138, t142, t143, t154, t157, t158, t159, t16, t160, t162, t165, &
811 : t167, t168, t17, t176, t180, t181, t184, t185, t188, t19, t190, t191, t195, t196, t199, &
812 : t2, t20, t202, t203, t204, t207, t208, t21, t210, t211, t22, t23, t237, t24, t240, t245, &
813 : t248, t249, t25, t255, t256, t258, t26, t265, t267, t272, t285, t288, t289, t29, t297, &
814 : t298, t3, t30, t301, t305, t309, t31, t313, t317, t32, t320, t321, t33, t338, t34, t341, &
815 : t35, t356, t36, t37, t376, t377, t382, t383, t387, t388, t391
816 : REAL(KIND=dp) :: t392, t395, t396, t399, t4, t400, t403, t404, t41, t416, t419, t42, t43, &
817 : t434, t458, t459, t47, t471, t472, t48, t484, t487, t49, t5, t50, t502, t51, t54, t57, &
818 : t58, t59, t6, t60, t62, t63, t66, t67, t68, t69, t70, t71, t72, t76, t77, t78, t79, t81, &
819 : t82, t85, t86, t87, t9, t90, t93, t96, t99, yval
820 :
821 2890201 : t1 = pi**(0.1e1_dp/0.3e1_dp)
822 2890201 : t2 = t1**2
823 2890201 : t3 = rho**(0.1e1_dp/0.3e1_dp)
824 2890201 : t4 = t3**2
825 2890201 : t5 = t4*rho
826 2890201 : t6 = t2*t5
827 2890201 : t9 = ndrho**2
828 2890201 : t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9/rho/0.4e1_dp)/0.3e1_dp
829 2890201 : t17 = 0.1e1_dp/t16
830 2890201 : yval = 0.2e1_dp/0.3e1_dp*t6*t17
831 2890201 : t19 = t3*rho
832 2890201 : t20 = 0.3141592654e1_dp**(0.1e1_dp/0.3e1_dp)
833 2890201 : t21 = t19*t20
834 2890201 : t22 = 0.1e1_dp/br_BB
835 2890201 : t23 = 0.1e1_dp/t2
836 2890201 : t24 = t22*t23
837 2890201 : t25 = 0.1e1_dp/t5
838 2890201 : t26 = t25*t16
839 2890201 : t29 = br_BB**2
840 2890201 : t30 = t1*pi
841 2890201 : t31 = t29*t30
842 2890201 : t32 = rho**2
843 2890201 : t33 = t32*rho
844 2890201 : t34 = t3*t33
845 2890201 : t35 = t16**2
846 2890201 : t36 = 0.1e1_dp/t35
847 2890201 : t37 = t34*t36
848 2890201 : t41 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t31*t37)
849 2890201 : t42 = t41*t22
850 2890201 : t43 = t23*t25
851 : t47 = 0.1500000000000000e1_dp*t24*t26 + 0.3e1_dp/0.2e1_dp*t42*t43 &
852 2890201 : *t16
853 2890201 : t48 = LOG(t47)
854 2890201 : t49 = t48 + 0.2e1_dp
855 2890201 : t50 = br_d1*t2
856 2890201 : t51 = t5*t17
857 2890201 : t54 = br_d2*t30
858 2890201 : t57 = pi**2
859 2890201 : t58 = br_d3*t57
860 2890201 : t59 = t32**2
861 2890201 : t60 = t59*rho
862 2890201 : t62 = 0.1e1_dp/t35/t16
863 2890201 : t63 = t60*t62
864 2890201 : t66 = t2*t57
865 2890201 : t67 = br_d4*t66
866 2890201 : t68 = t59*t32
867 2890201 : t69 = t4*t68
868 2890201 : t70 = t35**2
869 2890201 : t71 = 0.1e1_dp/t70
870 2890201 : t72 = t69*t71
871 2890201 : t76 = t1*t57*pi
872 2890201 : t77 = br_d5*t76
873 2890201 : t78 = t59**2
874 2890201 : t79 = t3*t78
875 2890201 : t81 = 0.1e1_dp/t70/t16
876 2890201 : t82 = t79*t81
877 : t85 = br_d0 + 0.2e1_dp/0.3e1_dp*t50*t51 + 0.4e1_dp/0.9e1_dp*t54*t37 &
878 : + 0.8e1_dp/0.27e2_dp*t58*t63 + 0.16e2_dp/0.81e2_dp*t67*t72 + 0.32e2_dp &
879 2890201 : /0.243e3_dp*t77*t82
880 2890201 : t86 = t49*t85
881 2890201 : t87 = br_e1*t2
882 2890201 : t90 = br_e2*t30
883 2890201 : t93 = br_e3*t57
884 2890201 : t96 = br_e4*t66
885 2890201 : t99 = br_e5*t76
886 : t102 = br_e0 + 0.2e1_dp/0.3e1_dp*t87*t51 + 0.4e1_dp/0.9e1_dp*t90*t37 &
887 : + 0.8e1_dp/0.27e2_dp*t93*t63 + 0.16e2_dp/0.81e2_dp*t96*t72 + 0.32e2_dp &
888 2890201 : /0.243e3_dp*t99*t82
889 2890201 : t103 = 0.1e1_dp/t102
890 2890201 : t104 = t86*t103
891 2890201 : t106 = EXP(t104/0.3e1_dp)
892 2890201 : t107 = t21*t106
893 2890201 : t108 = 0.1e1_dp/t49
894 2890201 : t109 = 0.1e1_dp/t85
895 2890201 : t110 = t108*t109
896 2890201 : t111 = EXP(-t104)
897 2890201 : t112 = t103*t111
898 2890201 : t115 = 0.1e1_dp - t111 - t86*t112/0.2e1_dp
899 2890201 : t117 = t110*t102*t115
900 2890201 : e_0 = e_0 + (-t107*t117)*sx
901 2890201 : IF (grad_deriv >= 1) THEN
902 2890201 : t124 = 0.1e1_dp/t4/t32
903 2890201 : t131 = 0.1e1_dp/t4/t33*gamma*t9
904 2890201 : t134 = 0.1e1_dp/t41
905 2890201 : t137 = t3*t32
906 2890201 : t138 = t137*t36
907 2890201 : t142 = t62*gamma
908 2890201 : t143 = t142*t9
909 2890201 : t154 = t42*t23
910 : t157 = -0.2500000000e1_dp*t24*t124*t16 - 0.1250000000e0_dp*t24* &
911 : t131 + 0.3e1_dp/0.4e1_dp*t134*t22*t23*t26*(0.40e2_dp/0.27e2_dp* &
912 : t31*t138 + 0.2e1_dp/0.27e2_dp*t31*t19*t143) - 0.5e1_dp/0.2e1_dp* &
913 2890201 : t42*t23*t124*t16 - t154*t131/0.8e1_dp
914 2890201 : t158 = 0.1e1_dp/t47
915 2890201 : t159 = t157*t158
916 2890201 : t160 = t85*t103
917 2890201 : t162 = t4*t17
918 2890201 : t165 = 0.1e1_dp/t3
919 2890201 : t167 = t36*gamma
920 2890201 : t168 = t167*t9
921 2890201 : t176 = t59*t62
922 2890201 : t180 = t71*gamma
923 2890201 : t181 = t180*t9
924 2890201 : t184 = t4*t60
925 2890201 : t185 = t184*t71
926 2890201 : t188 = t4*t59
927 2890201 : t190 = t81*gamma
928 2890201 : t191 = t190*t9
929 2890201 : t195 = t3*t59*t33
930 2890201 : t196 = t195*t81
931 2890201 : t199 = t3*t68
932 2890201 : t202 = 0.1e1_dp/t70/t35
933 2890201 : t203 = t202*gamma
934 2890201 : t204 = t203*t9
935 : t207 = 0.10e2_dp/0.9e1_dp*t50*t162 + t50*t165*t168/0.18e2_dp + 0.40e2_dp &
936 : /0.27e2_dp*t54*t138 + 0.2e1_dp/0.27e2_dp*t54*t19*t143 + &
937 : 0.40e2_dp/0.27e2_dp*t58*t176 + 0.2e1_dp/0.27e2_dp*t58*t33*t181 + &
938 : 0.320e3_dp/0.243e3_dp*t67*t185 + 0.16e2_dp/0.243e3_dp*t67*t188* &
939 : t191 + 0.800e3_dp/0.729e3_dp*t77*t196 + 0.40e2_dp/0.729e3_dp*t77*t199 &
940 2890201 : *t204
941 2890201 : t208 = t49*t207
942 2890201 : t210 = t102**2
943 2890201 : t211 = 0.1e1_dp/t210
944 : t237 = 0.10e2_dp/0.9e1_dp*t87*t162 + t87*t165*t168/0.18e2_dp + 0.40e2_dp &
945 : /0.27e2_dp*t90*t138 + 0.2e1_dp/0.27e2_dp*t90*t19*t143 + &
946 : 0.40e2_dp/0.27e2_dp*t93*t176 + 0.2e1_dp/0.27e2_dp*t93*t33*t181 + &
947 : 0.320e3_dp/0.243e3_dp*t96*t185 + 0.16e2_dp/0.243e3_dp*t96*t188* &
948 : t191 + 0.800e3_dp/0.729e3_dp*t99*t196 + 0.40e2_dp/0.729e3_dp*t99*t199 &
949 2890201 : *t204
950 2890201 : t240 = t159*t160 + t208*t103 - t86*t211*t237
951 2890201 : t245 = t49**2
952 2890201 : t248 = t21*t106/t245
953 2890201 : t249 = t109*t102
954 2890201 : t255 = t21*t106*t108
955 2890201 : t256 = t85**2
956 2890201 : t258 = 0.1e1_dp/t256*t102
957 2890201 : t265 = -t240
958 2890201 : t267 = t160*t111
959 2890201 : t272 = t211*t111
960 : e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t3*t20*t106*t117 - t21*t240*t106 &
961 : *t117/0.3e1_dp + t248*t249*t115*t157*t158 + t255*t258* &
962 : t115*t207 - t107*t110*t237*t115 - t107*t110*t102*(-t265 &
963 : *t111 - t159*t267/0.2e1_dp - t208*t112/0.2e1_dp + t86*t272 &
964 2890201 : *t237/0.2e1_dp - t86*t103*t265*t111/0.2e1_dp))*sx
965 2890201 : t285 = t124*gamma*ndrho
966 2890201 : t288 = t134*br_BB
967 2890201 : t289 = t288*t2
968 : t297 = 0.2500000000000000e0_dp*t24*t285 - t289*t4*t36*gamma* &
969 2890201 : ndrho/0.9e1_dp + t154*t285/0.4e1_dp
970 2890201 : t298 = t297*t158
971 2890201 : t301 = t167*ndrho
972 2890201 : t305 = t142*ndrho
973 2890201 : t309 = t180*ndrho
974 2890201 : t313 = t190*ndrho
975 2890201 : t317 = t203*ndrho
976 : t320 = -t50*t4*t301/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t54*t137*t305 &
977 : - 0.4e1_dp/0.27e2_dp*t58*t59*t309 - 0.32e2_dp/0.243e3_dp*t67*t184 &
978 2890201 : *t313 - 0.80e2_dp/0.729e3_dp*t77*t195*t317
979 2890201 : t321 = t49*t320
980 : t338 = -t87*t4*t301/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t90*t137*t305 &
981 : - 0.4e1_dp/0.27e2_dp*t93*t59*t309 - 0.32e2_dp/0.243e3_dp*t96*t184 &
982 2890201 : *t313 - 0.80e2_dp/0.729e3_dp*t99*t195*t317
983 2890201 : t341 = t298*t160 + t321*t103 - t86*t211*t338
984 2890201 : t356 = -t341
985 : e_ndrho = e_ndrho + (-t21*t341*t106*t117/0.3e1_dp + t248*t249*t115*t297 &
986 : *t158 + t255*t258*t115*t320 - t107*t110*t338*t115 &
987 : - t107*t110*t102*(-t356*t111 - t298*t267/0.2e1_dp - t321 &
988 : *t112/0.2e1_dp + t86*t272*t338/0.2e1_dp - t86*t103*t356* &
989 2890201 : t111/0.2e1_dp))*sx
990 2890201 : t376 = t5*t36
991 2890201 : t377 = t376*gamma
992 : t382 = -0.1000000000e1_dp*t24*t25*gamma + 0.4e1_dp/0.9e1_dp*t289*t377 &
993 2890201 : - t42*t43*gamma
994 2890201 : t383 = t382*t158
995 2890201 : t387 = t34*t62
996 2890201 : t388 = t387*gamma
997 2890201 : t391 = t60*t71
998 2890201 : t392 = t391*gamma
999 2890201 : t395 = t69*t81
1000 2890201 : t396 = t395*gamma
1001 2890201 : t399 = t79*t202
1002 2890201 : t400 = t399*gamma
1003 : t403 = 0.4e1_dp/0.9e1_dp*t50*t377 + 0.16e2_dp/0.27e2_dp*t54*t388 + &
1004 : 0.16e2_dp/0.27e2_dp*t58*t392 + 0.128e3_dp/0.243e3_dp*t67*t396 + 0.320e3_dp &
1005 2890201 : /0.729e3_dp*t77*t400
1006 2890201 : t404 = t49*t403
1007 : t416 = 0.4e1_dp/0.9e1_dp*t87*t377 + 0.16e2_dp/0.27e2_dp*t90*t388 + &
1008 : 0.16e2_dp/0.27e2_dp*t93*t392 + 0.128e3_dp/0.243e3_dp*t96*t396 + 0.320e3_dp &
1009 2890201 : /0.729e3_dp*t99*t400
1010 2890201 : t419 = t383*t160 + t404*t103 - t86*t211*t416
1011 2890201 : t434 = -t419
1012 : e_tau = e_tau + (-t21*t419*t106*t117/0.3e1_dp + t248*t249*t115*t382 &
1013 : *t158 + t255*t258*t115*t403 - t107*t110*t416*t115 - &
1014 : t107*t110*t102*(-t434*t111 - t383*t267/0.2e1_dp - t404* &
1015 : t112/0.2e1_dp + t86*t272*t416/0.2e1_dp - t86*t103*t434*t111 &
1016 2890201 : /0.2e1_dp))*sx
1017 : t458 = 0.2500000000000000e0_dp*t24*t25 - t288*t6*t36/0.9e1_dp + &
1018 2890201 : t42*t43/0.4e1_dp
1019 2890201 : t459 = t458*t158
1020 : t471 = -t50*t376/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t54*t387 - 0.4e1_dp/ &
1021 : 0.27e2_dp*t58*t391 - 0.32e2_dp/0.243e3_dp*t67*t395 - 0.80e2_dp/0.729e3_dp &
1022 2890201 : *t77*t399
1023 2890201 : t472 = t49*t471
1024 : t484 = -t87*t376/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t90*t387 - 0.4e1_dp/ &
1025 : 0.27e2_dp*t93*t391 - 0.32e2_dp/0.243e3_dp*t96*t395 - 0.80e2_dp/0.729e3_dp &
1026 2890201 : *t99*t399
1027 2890201 : t487 = t459*t160 + t472*t103 - t86*t211*t484
1028 2890201 : t502 = -t487
1029 : e_laplace_rho = e_laplace_rho + (-t21*t487*t106*t117/0.3e1_dp + t248*t249*t115*t458 &
1030 : *t158 + t255*t258*t115*t471 - t107*t110*t484*t115 - &
1031 : t107*t110*t102*(-t502*t111 - t459*t267/0.2e1_dp - t472* &
1032 : t112/0.2e1_dp + t86*t272*t484/0.2e1_dp - t86*t103*t502*t111 &
1033 2890201 : /0.2e1_dp))*sx
1034 : END IF
1035 :
1036 2890201 : END SUBROUTINE x_br_lsd_y_gt_0
1037 :
1038 : ! **************************************************************************************************
1039 : !> \brief Truncated long range part for y <= 0
1040 : !> \param rho ...
1041 : !> \param ndrho ...
1042 : !> \param tau ...
1043 : !> \param laplace_rho ...
1044 : !> \param e_0 ...
1045 : !> \param e_rho ...
1046 : !> \param e_ndrho ...
1047 : !> \param e_tau ...
1048 : !> \param e_laplace_rho ...
1049 : !> \param sx ...
1050 : !> \param R ...
1051 : !> \param gamma ...
1052 : !> \param grad_deriv ...
1053 : !> \par History
1054 : !> 12.2008 created [mguidon]
1055 : !> \author mguidon
1056 : ! **************************************************************************************************
1057 386143 : SUBROUTINE x_br_lsd_y_lte_0_cutoff(rho, ndrho, tau, laplace_rho, e_0, &
1058 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1059 : sx, R, gamma, grad_deriv)
1060 : REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho
1061 : REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho
1062 : REAL(dp), INTENT(IN) :: sx, R, gamma
1063 : INTEGER, INTENT(IN) :: grad_deriv
1064 :
1065 : REAL(KIND=dp) :: brval, t1, t101, t11, t12, t18, t2, t20, &
1066 : t24, t25, t26, t3, t31, t33, t36, t38, &
1067 : t4, t41, t43, t47, t50, t54, t56, t6, &
1068 : t60, t62, t66, t69, t7, t70, t88, t89, &
1069 : t96
1070 :
1071 386143 : t1 = 8**(0.1e1_dp/0.3e1_dp)
1072 386143 : t2 = t1**2
1073 386143 : t3 = pi**(0.1e1_dp/0.3e1_dp)
1074 386143 : t4 = t3**2
1075 386143 : t6 = rho**(0.1e1_dp/0.3e1_dp)
1076 386143 : t7 = t6**2
1077 386143 : t11 = ndrho**2
1078 386143 : t12 = 0.1e1_dp/rho
1079 386143 : t18 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t11*t12/0.4e1_dp)/0.3e1_dp
1080 386143 : t20 = t7*rho/t18
1081 386143 : t24 = ATAN(0.2e1_dp/0.3e1_dp*br_a1*t4*t20 + br_a2)
1082 386143 : t25 = -t24 + br_a3
1083 386143 : t26 = t25**2
1084 386143 : t31 = t3*pi
1085 386143 : t33 = rho**2
1086 386143 : t36 = t18**2
1087 386143 : t38 = t6*t33*rho/t36
1088 386143 : t41 = pi**2
1089 386143 : t43 = t33**2
1090 386143 : t47 = t43*rho/t36/t18
1091 386143 : t50 = t4*t41
1092 386143 : t54 = t36**2
1093 386143 : t56 = t7*t43*t33/t54
1094 386143 : t60 = t3*t41*pi
1095 386143 : t62 = t43**2
1096 386143 : t66 = t6*t62/t54/t18
1097 : t69 = br_c0 + 0.2e1_dp/0.3e1_dp*br_c1*t4*t20 + 0.4e1_dp/0.9e1_dp*br_c2 &
1098 : *t31*t38 + 0.8e1_dp/0.27e2_dp*br_c3*t41*t47 + 0.16e2_dp/0.81e2_dp &
1099 386143 : *br_c4*t50*t56 + 0.32e2_dp/0.243e3_dp*br_c5*t60*t66
1100 386143 : t70 = t69**2
1101 : t88 = br_b0 + 0.2e1_dp/0.3e1_dp*br_b1*t4*t20 + 0.4e1_dp/0.9e1_dp*br_b2 &
1102 : *t31*t38 + 0.8e1_dp/0.27e2_dp*br_b3*t41*t47 + 0.16e2_dp/0.81e2_dp &
1103 386143 : *br_b4*t50*t56 + 0.32e2_dp/0.243e3_dp*br_b5*t60*t66
1104 386143 : t89 = t88**2
1105 386143 : t96 = EXP(-t25*t69/t88)
1106 : t101 = (t26*t25*t70*t69/t89/t88*t96/0.3141592654e1_dp* &
1107 386143 : t12)**(0.1e1_dp/0.3e1_dp)
1108 386143 : brval = REAL(t2, KIND=dp)*t101/0.8e1_dp
1109 :
1110 386143 : IF (R > brval) THEN
1111 : CALL x_br_lsd_y_lte_0_cutoff_R_gt_b(rho, ndrho, tau, laplace_rho, e_0, &
1112 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1113 368598 : sx, R, gamma, grad_deriv)
1114 : ELSE
1115 : CALL x_br_lsd_y_lte_0_cutoff_R_lte_b(rho, ndrho, tau, laplace_rho, e_0, &
1116 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1117 17545 : sx, R, gamma, grad_deriv)
1118 : END IF
1119 :
1120 386143 : END SUBROUTINE x_br_lsd_y_lte_0_cutoff
1121 :
1122 : ! **************************************************************************************************
1123 : !> \brief Truncated long range part for y <= 0 and R > b
1124 : !> \param rho ...
1125 : !> \param ndrho ...
1126 : !> \param tau ...
1127 : !> \param laplace_rho ...
1128 : !> \param e_0 ...
1129 : !> \param e_rho ...
1130 : !> \param e_ndrho ...
1131 : !> \param e_tau ...
1132 : !> \param e_laplace_rho ...
1133 : !> \param sx ...
1134 : !> \param R ...
1135 : !> \param gamma ...
1136 : !> \param grad_deriv ...
1137 : !> \par History
1138 : !> 12.2008 created [mguidon]
1139 : !> \author mguidon
1140 : ! **************************************************************************************************
1141 368598 : SUBROUTINE x_br_lsd_y_lte_0_cutoff_R_gt_b(rho, ndrho, tau, laplace_rho, e_0, &
1142 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1143 : sx, R, gamma, grad_deriv)
1144 : REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho
1145 : REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho
1146 : REAL(dp), INTENT(IN) :: sx, R, gamma
1147 : INTEGER, INTENT(IN) :: grad_deriv
1148 :
1149 : REAL(KIND=dp) :: t1, t10, t100, t101, t102, t103, t104, t106, t108, t109, t110, t114, t117, &
1150 : t119, t120, t123, t124, t129, t132, t134, t135, t138, t139, t141, t142, t143, t149, t150, &
1151 : t153, t155, t156, t159, t16, t163, t164, t167, t168, t17, t171, t173, t174, t178, t179, &
1152 : t18, t182, t185, t186, t187, t190, t192, t193, t2, t21, t219, t22, t221, t223, t224, &
1153 : t225, t227, t228, t23, t231, t232, t233, t235, t24, t240, t245, t247, t259, t261, t263, &
1154 : t265, t268, t269, t27, t270, t272, t274, t275, t28, t283, t288, t29, t291, t3, t30, t301, &
1155 : t305, t31, t312, t314, t315, t317, t319, t32, t323, t327, t33
1156 : REAL(KIND=dp) :: t331, t335, t338, t34, t340, t356, t358, t359, t361, t363, t364, t366, &
1157 : t367, t37, t38, t388, t39, t390, t392, t394, t397, t399, t4, t40, t403, t405, t409, t410, &
1158 : t414, t42, t423, t426, t43, t434, t443, t450, t455, t456, t459, t46, t460, t463, t464, &
1159 : t467, t468, t47, t471, t472, t475, t477, t48, t488, t49, t490, t493, t494, t496, t497, &
1160 : t499, t5, t50, t51, t516, t518, t52, t520, t522, t525, t526, t527, t530, t532, t545, &
1161 : t548, t56, t563, t57, t572, t574, t58, t585, t587, t59, t598, t6, t600, t605, t606, t608, &
1162 : t609, t61, t62, t628, t630, t632, t634, t639, t641, t645, t65, t655
1163 : REAL(KIND=dp) :: t66, t67, t672, t70, t73, t76, t79, t82, t83, t84, t85, t86, t87, t88, t89, &
1164 : t9, t90, t91, t93, t94, t95, t96, t97, t99
1165 :
1166 368598 : t1 = pi**(0.1e1_dp/0.3e1_dp)
1167 368598 : t2 = t1**2
1168 368598 : t3 = br_a1*t2
1169 368598 : t4 = rho**(0.1e1_dp/0.3e1_dp)
1170 368598 : t5 = t4**2
1171 368598 : t6 = t5*rho
1172 368598 : t9 = ndrho**2
1173 368598 : t10 = 0.1e1_dp/rho
1174 368598 : t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9*t10/0.4e1_dp)/0.3e1_dp
1175 368598 : t17 = 0.1e1_dp/t16
1176 368598 : t18 = t6*t17
1177 368598 : t21 = 0.2e1_dp/0.3e1_dp*t3*t18 + br_a2
1178 368598 : t22 = ATAN(t21)
1179 368598 : t23 = -t22 + br_a3
1180 368598 : t24 = br_c1*t2
1181 368598 : t27 = t1*pi
1182 368598 : t28 = br_c2*t27
1183 368598 : t29 = rho**2
1184 368598 : t30 = t29*rho
1185 368598 : t31 = t4*t30
1186 368598 : t32 = t16**2
1187 368598 : t33 = 0.1e1_dp/t32
1188 368598 : t34 = t31*t33
1189 368598 : t37 = pi**2
1190 368598 : t38 = br_c3*t37
1191 368598 : t39 = t29**2
1192 368598 : t40 = t39*rho
1193 368598 : t42 = 0.1e1_dp/t32/t16
1194 368598 : t43 = t40*t42
1195 368598 : t46 = t2*t37
1196 368598 : t47 = br_c4*t46
1197 368598 : t48 = t39*t29
1198 368598 : t49 = t5*t48
1199 368598 : t50 = t32**2
1200 368598 : t51 = 0.1e1_dp/t50
1201 368598 : t52 = t49*t51
1202 368598 : t56 = t1*t37*pi
1203 368598 : t57 = br_c5*t56
1204 368598 : t58 = t39**2
1205 368598 : t59 = t4*t58
1206 368598 : t61 = 0.1e1_dp/t50/t16
1207 368598 : t62 = t59*t61
1208 : t65 = br_c0 + 0.2e1_dp/0.3e1_dp*t24*t18 + 0.4e1_dp/0.9e1_dp*t28*t34 &
1209 : + 0.8e1_dp/0.27e2_dp*t38*t43 + 0.16e2_dp/0.81e2_dp*t47*t52 + 0.32e2_dp &
1210 368598 : /0.243e3_dp*t57*t62
1211 368598 : t66 = t23*t65
1212 368598 : t67 = br_b1*t2
1213 368598 : t70 = br_b2*t27
1214 368598 : t73 = br_b3*t37
1215 368598 : t76 = br_b4*t46
1216 368598 : t79 = br_b5*t56
1217 : t82 = br_b0 + 0.2e1_dp/0.3e1_dp*t67*t18 + 0.4e1_dp/0.9e1_dp*t70*t34 &
1218 : + 0.8e1_dp/0.27e2_dp*t73*t43 + 0.16e2_dp/0.81e2_dp*t76*t52 + 0.32e2_dp &
1219 368598 : /0.243e3_dp*t79*t62
1220 368598 : t83 = 0.1e1_dp/t82
1221 368598 : t84 = t66*t83
1222 368598 : t85 = 8**(0.1e1_dp/0.3e1_dp)
1223 368598 : t86 = t23**2
1224 368598 : t87 = t86*t23
1225 368598 : t88 = t65**2
1226 368598 : t89 = t88*t65
1227 368598 : t90 = t87*t89
1228 368598 : t91 = t82**2
1229 368598 : t93 = 0.1e1_dp/t91/t82
1230 368598 : t94 = t90*t93
1231 368598 : t95 = EXP(-t84)
1232 368598 : t96 = 0.1e1_dp/0.3141592654e1_dp
1233 368598 : t97 = t95*t96
1234 368598 : t99 = t94*t97*t10
1235 368598 : t100 = t99**(0.1e1_dp/0.3e1_dp)
1236 368598 : t101 = 0.1e1_dp/t100
1237 368598 : t102 = REAL(t85, KIND=dp)*t101
1238 368598 : t103 = t102*R
1239 368598 : t104 = t84*t103
1240 368598 : t106 = EXP(t84 - t104)
1241 368598 : t108 = t106*t23
1242 368598 : t109 = t65*t83
1243 368598 : t110 = t108*t109
1244 368598 : t114 = t83*REAL(t85, KIND=dp)*t101*R
1245 368598 : t117 = EXP(-t84 - t104)
1246 368598 : t119 = t117*t23
1247 368598 : t120 = t119*t109
1248 : t123 = -0.2e1_dp*t106 + t110 - t108*t65*t114 + 0.2e1_dp*t117 + t120 &
1249 368598 : + t119*t65*t114
1250 368598 : t124 = rho*t123
1251 368598 : e_0 = e_0 + (t124*t102/0.8e1_dp)*sx
1252 368598 : IF (grad_deriv >= 1) THEN
1253 194202 : t129 = t5*t17
1254 194202 : t132 = 0.1e1_dp/t4
1255 194202 : t134 = t33*gamma
1256 194202 : t135 = t134*t9
1257 194202 : t138 = 0.10e2_dp/0.9e1_dp*t3*t129 + t3*t132*t135/0.18e2_dp
1258 194202 : t139 = t21**2
1259 194202 : t141 = 0.1e1_dp/(0.1e1_dp + t139)
1260 194202 : t142 = t138*t141
1261 194202 : t143 = t142*t109
1262 194202 : t149 = t4*t29
1263 194202 : t150 = t149*t33
1264 194202 : t153 = t4*rho
1265 194202 : t155 = t42*gamma
1266 194202 : t156 = t155*t9
1267 194202 : t159 = t39*t42
1268 194202 : t163 = t51*gamma
1269 194202 : t164 = t163*t9
1270 194202 : t167 = t5*t40
1271 194202 : t168 = t167*t51
1272 194202 : t171 = t5*t39
1273 194202 : t173 = t61*gamma
1274 194202 : t174 = t173*t9
1275 194202 : t178 = t4*t39*t30
1276 194202 : t179 = t178*t61
1277 194202 : t182 = t4*t48
1278 194202 : t185 = 0.1e1_dp/t50/t32
1279 194202 : t186 = t185*gamma
1280 194202 : t187 = t186*t9
1281 : t190 = 0.10e2_dp/0.9e1_dp*t24*t129 + t24*t132*t135/0.18e2_dp + 0.40e2_dp &
1282 : /0.27e2_dp*t28*t150 + 0.2e1_dp/0.27e2_dp*t28*t153*t156 + &
1283 : 0.40e2_dp/0.27e2_dp*t38*t159 + 0.2e1_dp/0.27e2_dp*t38*t30*t164 &
1284 : + 0.320e3_dp/0.243e3_dp*t47*t168 + 0.16e2_dp/0.243e3_dp*t47*t171* &
1285 : t174 + 0.800e3_dp/0.729e3_dp*t57*t179 + 0.40e2_dp/0.729e3_dp*t57* &
1286 194202 : t182*t187
1287 194202 : t192 = t23*t190*t83
1288 194202 : t193 = 0.1e1_dp/t91
1289 : t219 = 0.10e2_dp/0.9e1_dp*t67*t129 + t67*t132*t135/0.18e2_dp + 0.40e2_dp &
1290 : /0.27e2_dp*t70*t150 + 0.2e1_dp/0.27e2_dp*t70*t153*t156 + &
1291 : 0.40e2_dp/0.27e2_dp*t73*t159 + 0.2e1_dp/0.27e2_dp*t73*t30*t164 &
1292 : + 0.320e3_dp/0.243e3_dp*t76*t168 + 0.16e2_dp/0.243e3_dp*t76*t171* &
1293 : t174 + 0.800e3_dp/0.729e3_dp*t79*t179 + 0.40e2_dp/0.729e3_dp*t79* &
1294 194202 : t182*t187
1295 194202 : t221 = t66*t193*t219
1296 194202 : t223 = t142*t65*t114
1297 194202 : t224 = t192*t103
1298 194202 : t225 = t66*t193
1299 194202 : t227 = t102*R*t219
1300 194202 : t228 = t225*t227
1301 194202 : t231 = REAL(t85, KIND=dp)/t100/t99
1302 194202 : t232 = t86*t89
1303 194202 : t233 = t93*t95
1304 194202 : t235 = t96*t10
1305 194202 : t240 = t87*t88*t93
1306 194202 : t245 = t91**2
1307 194202 : t247 = t90/t245
1308 : t259 = -0.3e1_dp*t232*t233*t235*t142 + 0.3e1_dp*t240*t97*t10 &
1309 : *t190 - 0.3e1_dp*t247*t97*t10*t219 + t94*(t143 - t192 + &
1310 194202 : t221)*t95*t235 - t94*t97/t29
1311 194202 : t261 = t231*R*t259
1312 194202 : t263 = t84*t261/0.3e1_dp
1313 194202 : t265 = (-t143 + t192 - t221 + t223 - t224 + t228 + t263)*t106
1314 194202 : t268 = t106*t138
1315 194202 : t269 = t141*t65
1316 194202 : t270 = t269*t83
1317 194202 : t272 = t190*t83
1318 194202 : t274 = t65*t193
1319 194202 : t275 = t274*t219
1320 194202 : t283 = t108*t274
1321 194202 : t288 = (t143 - t192 + t221 + t223 - t224 + t228 + t263)*t117
1322 194202 : t291 = t117*t138
1323 194202 : t301 = t119*t274
1324 : t305 = -0.2e1_dp*t265 + t265*t84 - t268*t270 + t108*t272 - t108 &
1325 : *t275 - t265*t66*t114 + t268*t269*t114 - t108*t190* &
1326 : t114 + t283*t227 + t110*t261/0.3e1_dp + 0.2e1_dp*t288 + t288*t84 &
1327 : - t291*t270 + t119*t272 - t119*t275 + t288*t66*t114 - &
1328 : t291*t269*t114 + t119*t190*t114 - t301*t227 - t120*t261 &
1329 194202 : /0.3e1_dp
1330 : e_rho = e_rho + (t123*REAL(t85, KIND=dp)*t101/0.8e1_dp + rho*t305*t102/0.8e1_dp &
1331 194202 : - t124*t231*t259/0.24e2_dp)*sx
1332 194202 : t312 = t5*t33
1333 194202 : t314 = gamma*ndrho
1334 194202 : t315 = t314*t270
1335 194202 : t317 = t3*t312*t315/0.9e1_dp
1336 194202 : t319 = t134*ndrho
1337 194202 : t323 = t155*ndrho
1338 194202 : t327 = t163*ndrho
1339 194202 : t331 = t173*ndrho
1340 194202 : t335 = t186*ndrho
1341 : t338 = -t24*t5*t319/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t28*t149*t323 &
1342 : - 0.4e1_dp/0.27e2_dp*t38*t39*t327 - 0.32e2_dp/0.243e3_dp*t47*t167 &
1343 194202 : *t331 - 0.80e2_dp/0.729e3_dp*t57*t178*t335
1344 194202 : t340 = t23*t338*t83
1345 : t356 = -t67*t5*t319/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t70*t149*t323 &
1346 : - 0.4e1_dp/0.27e2_dp*t73*t39*t327 - 0.32e2_dp/0.243e3_dp*t76*t167 &
1347 194202 : *t331 - 0.80e2_dp/0.729e3_dp*t79*t178*t335
1348 194202 : t358 = t66*t193*t356
1349 194202 : t359 = t3*t5
1350 194202 : t361 = t270*t103
1351 194202 : t363 = t359*t319*t361/0.9e1_dp
1352 194202 : t364 = t340*t103
1353 194202 : t366 = t102*R*t356
1354 194202 : t367 = t225*t366
1355 : t388 = t232*t93*t97*t132*t3*t33*t314*t141/0.3e1_dp + 0.3e1_dp &
1356 : *t240*t97*t10*t338 - 0.3e1_dp*t247*t97*t10*t356 + &
1357 194202 : t94*(-t317 - t340 + t358)*t95*t235
1358 194202 : t390 = t231*R*t388
1359 194202 : t392 = t84*t390/0.3e1_dp
1360 194202 : t394 = (t317 + t340 - t358 - t363 - t364 + t367 + t392)*t106
1361 194202 : t397 = t106*br_a1
1362 194202 : t399 = t2*t5*t33
1363 194202 : t403 = t338*t83
1364 194202 : t405 = t274*t356
1365 194202 : t409 = t397*t2
1366 194202 : t410 = t312*gamma
1367 194202 : t414 = ndrho*t141*t65*t114
1368 194202 : t423 = (-t317 - t340 + t358 - t363 - t364 + t367 + t392)*t117
1369 194202 : t426 = t117*br_a1
1370 194202 : t434 = t426*t2
1371 : t443 = -0.2e1_dp*t394 + t394*t84 + t397*t399*t315/0.9e1_dp + t108 &
1372 : *t403 - t108*t405 - t394*t66*t114 - t409*t410*t414/ &
1373 : 0.9e1_dp - t108*t338*t114 + t283*t366 + t110*t390/0.3e1_dp + &
1374 : 0.2e1_dp*t423 + t423*t84 + t426*t399*t315/0.9e1_dp + t119*t403 &
1375 : - t119*t405 + t423*t66*t114 + t434*t410*t414/0.9e1_dp &
1376 194202 : + t119*t338*t114 - t301*t366 - t120*t390/0.3e1_dp
1377 194202 : e_ndrho = e_ndrho + (rho*t443*t102/0.8e1_dp - t124*t231*t388/0.24e2_dp)*sx
1378 194202 : t450 = t6*t33
1379 194202 : t455 = 0.4e1_dp/0.9e1_dp*t3*t450*gamma*t141*t109
1380 194202 : t456 = t450*gamma
1381 194202 : t459 = t31*t42
1382 194202 : t460 = t459*gamma
1383 194202 : t463 = t40*t51
1384 194202 : t464 = t463*gamma
1385 194202 : t467 = t49*t61
1386 194202 : t468 = t467*gamma
1387 194202 : t471 = t59*t185
1388 194202 : t472 = t471*gamma
1389 : t475 = 0.4e1_dp/0.9e1_dp*t24*t456 + 0.16e2_dp/0.27e2_dp*t28*t460 + &
1390 : 0.16e2_dp/0.27e2_dp*t38*t464 + 0.128e3_dp/0.243e3_dp*t47*t468 + 0.320e3_dp &
1391 194202 : /0.729e3_dp*t57*t472
1392 194202 : t477 = t23*t475*t83
1393 : t488 = 0.4e1_dp/0.9e1_dp*t67*t456 + 0.16e2_dp/0.27e2_dp*t70*t460 + &
1394 : 0.16e2_dp/0.27e2_dp*t73*t464 + 0.128e3_dp/0.243e3_dp*t76*t468 + 0.320e3_dp &
1395 194202 : /0.729e3_dp*t79*t472
1396 194202 : t490 = t66*t193*t488
1397 194202 : t493 = 0.4e1_dp/0.9e1_dp*t3*t456*t361
1398 194202 : t494 = t477*t103
1399 194202 : t496 = t102*R*t488
1400 194202 : t497 = t225*t496
1401 194202 : t499 = t232*t233*t96
1402 : t516 = -0.4e1_dp/0.3e1_dp*t499*t359*t134*t141 + 0.3e1_dp*t240* &
1403 : t97*t10*t475 - 0.3e1_dp*t247*t97*t10*t488 + t94*(t455 - &
1404 194202 : t477 + t490)*t95*t235
1405 194202 : t518 = t231*R*t516
1406 194202 : t520 = t84*t518/0.3e1_dp
1407 194202 : t522 = (-t455 + t477 - t490 + t493 - t494 + t497 + t520)*t106
1408 194202 : t525 = t2*t6
1409 194202 : t526 = t397*t525
1410 194202 : t527 = t134*t270
1411 194202 : t530 = t475*t83
1412 194202 : t532 = t274*t488
1413 194202 : t545 = (t455 - t477 + t490 + t493 - t494 + t497 + t520)*t117
1414 194202 : t548 = t426*t525
1415 : t563 = -0.2e1_dp*t522 + t522*t84 - 0.4e1_dp/0.9e1_dp*t526*t527 + t108 &
1416 : *t530 - t108*t532 - t522*t66*t114 + 0.4e1_dp/0.9e1_dp*t409 &
1417 : *t456*t361 - t108*t475*t114 + t283*t496 + t110*t518/ &
1418 : 0.3e1_dp + 0.2e1_dp*t545 + t545*t84 - 0.4e1_dp/0.9e1_dp*t548*t527 + &
1419 : t119*t530 - t119*t532 + t545*t66*t114 - 0.4e1_dp/0.9e1_dp*t434 &
1420 : *t456*t361 + t119*t475*t114 - t301*t496 - t120*t518 &
1421 194202 : /0.3e1_dp
1422 194202 : e_tau = e_tau + (rho*t563*t102/0.8e1_dp - t124*t231*t516/0.24e2_dp)*sx
1423 194202 : t572 = t33*t141*t109
1424 194202 : t574 = t3*t6*t572/0.9e1_dp
1425 : t585 = -t24*t450/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t28*t459 - 0.4e1_dp/ &
1426 : 0.27e2_dp*t38*t463 - 0.32e2_dp/0.243e3_dp*t47*t467 - 0.80e2_dp/0.729e3_dp &
1427 194202 : *t57*t471
1428 194202 : t587 = t23*t585*t83
1429 : t598 = -t67*t450/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t70*t459 - 0.4e1_dp/ &
1430 : 0.27e2_dp*t73*t463 - 0.32e2_dp/0.243e3_dp*t76*t467 - 0.80e2_dp/0.729e3_dp &
1431 194202 : *t79*t471
1432 194202 : t600 = t66*t193*t598
1433 194202 : t605 = t3*t450*t141*t109*t103/0.9e1_dp
1434 194202 : t606 = t587*t103
1435 194202 : t608 = t102*R*t598
1436 194202 : t609 = t225*t608
1437 : t628 = t499*t5*br_a1*t2*t33*t141/0.3e1_dp + 0.3e1_dp*t240* &
1438 : t97*t10*t585 - 0.3e1_dp*t247*t97*t10*t598 + t94*(-t574 &
1439 194202 : - t587 + t600)*t95*t235
1440 194202 : t630 = t231*R*t628
1441 194202 : t632 = t84*t630/0.3e1_dp
1442 194202 : t634 = (t574 + t587 - t600 - t605 - t606 + t609 + t632)*t106
1443 194202 : t639 = t585*t83
1444 194202 : t641 = t274*t598
1445 194202 : t645 = t525*t33
1446 194202 : t655 = (-t574 - t587 + t600 - t605 - t606 + t609 + t632)*t117
1447 : t672 = -0.2e1_dp*t634 + t634*t84 + t526*t572/0.9e1_dp + t108*t639 &
1448 : - t108*t641 - t634*t66*t114 - t397*t645*t361/0.9e1_dp &
1449 : - t108*t585*t114 + t283*t608 + t110*t630/0.3e1_dp + 0.2e1_dp* &
1450 : t655 + t655*t84 + t548*t572/0.9e1_dp + t119*t639 - t119*t641 &
1451 : + t655*t66*t114 + t426*t645*t361/0.9e1_dp + t119*t585 &
1452 194202 : *t114 - t301*t608 - t120*t630/0.3e1_dp
1453 194202 : e_laplace_rho = e_laplace_rho + (rho*t672*t102/0.8e1_dp - t124*t231*t628/0.24e2_dp)*sx
1454 : END IF
1455 :
1456 368598 : END SUBROUTINE x_br_lsd_y_lte_0_cutoff_R_gt_b
1457 :
1458 : ! **************************************************************************************************
1459 : !> \brief Truncated long range part for y <= 0 and R <= b
1460 : !> \param rho ...
1461 : !> \param ndrho ...
1462 : !> \param tau ...
1463 : !> \param laplace_rho ...
1464 : !> \param e_0 ...
1465 : !> \param e_rho ...
1466 : !> \param e_ndrho ...
1467 : !> \param e_tau ...
1468 : !> \param e_laplace_rho ...
1469 : !> \param sx ...
1470 : !> \param R ...
1471 : !> \param gamma ...
1472 : !> \param grad_deriv ...
1473 : !> \par History
1474 : !> 12.2008 created [mguidon]
1475 : !> \author mguidon
1476 : ! **************************************************************************************************
1477 17545 : SUBROUTINE x_br_lsd_y_lte_0_cutoff_R_lte_b(rho, ndrho, tau, laplace_rho, e_0, &
1478 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1479 : sx, R, gamma, grad_deriv)
1480 : REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho
1481 : REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho
1482 : REAL(dp), INTENT(IN) :: sx, R, gamma
1483 : INTEGER, INTENT(IN) :: grad_deriv
1484 :
1485 : REAL(KIND=dp) :: t1, t10, t100, t101, t102, t103, t104, t106, t108, t109, t110, t111, t113, &
1486 : t115, t116, t118, t119, t121, t123, t128, t131, t133, t134, t137, t138, t140, t141, t143, &
1487 : t146, t153, t154, t157, t159, t16, t160, t163, t167, t168, t17, t171, t172, t175, t177, &
1488 : t178, t18, t182, t183, t186, t189, t190, t191, t194, t196, t197, t199, t2, t200, t21, &
1489 : t22, t226, t229, t23, t232, t233, t234, t235, t237, t24, t242, t247, t249, t254, t256, &
1490 : t264, t267, t27, t270, t272, t277, t28, t280, t281, t29, t292, t295, t298, t299, t3, t30, &
1491 : t302, t31, t310, t314, t315, t317, t318, t32, t324, t328, t33
1492 : REAL(KIND=dp) :: t332, t336, t339, t34, t341, t342, t359, t362, t368, t369, t37, t38, t383, &
1493 : t385, t387, t39, t392, t395, t398, t4, t40, t402, t404, t409, t42, t420, t427, t428, &
1494 : t429, t43, t432, t443, t444, t446, t450, t451, t454, t455, t458, t459, t46, t462, t463, &
1495 : t466, t468, t469, t47, t48, t481, t484, t487, t49, t5, t50, t501, t504, t506, t51, t511, &
1496 : t514, t517, t52, t521, t525, t529, t540, t547, t548, t549, t552, t56, t566, t57, t578, &
1497 : t58, t580, t581, t59, t593, t596, t6, t61, t612, t613, t614, t616, t618, t62, t623, t626, &
1498 : t629, t638, t65, t654, t655, t656, t659, t66, t67, t70, t73, t76
1499 : REAL(KIND=dp) :: t79, t82, t83, t84, t85, t86, t87, t88, t89, t9, t90, t91, t93, t94, t95, &
1500 : t96, t97, t99
1501 :
1502 17545 : t1 = pi**(0.1e1_dp/0.3e1_dp)
1503 17545 : t2 = t1**2
1504 17545 : t3 = br_a1*t2
1505 17545 : t4 = rho**(0.1e1_dp/0.3e1_dp)
1506 17545 : t5 = t4**2
1507 17545 : t6 = t5*rho
1508 17545 : t9 = ndrho**2
1509 17545 : t10 = 0.1e1_dp/rho
1510 17545 : t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9*t10/0.4e1_dp)/0.3e1_dp
1511 17545 : t17 = 0.1e1_dp/t16
1512 17545 : t18 = t6*t17
1513 17545 : t21 = 0.2e1_dp/0.3e1_dp*t3*t18 + br_a2
1514 17545 : t22 = ATAN(t21)
1515 17545 : t23 = -t22 + br_a3
1516 17545 : t24 = br_c1*t2
1517 17545 : t27 = t1*pi
1518 17545 : t28 = br_c2*t27
1519 17545 : t29 = rho**2
1520 17545 : t30 = t29*rho
1521 17545 : t31 = t4*t30
1522 17545 : t32 = t16**2
1523 17545 : t33 = 0.1e1_dp/t32
1524 17545 : t34 = t31*t33
1525 17545 : t37 = pi**2
1526 17545 : t38 = br_c3*t37
1527 17545 : t39 = t29**2
1528 17545 : t40 = t39*rho
1529 17545 : t42 = 0.1e1_dp/t32/t16
1530 17545 : t43 = t40*t42
1531 17545 : t46 = t2*t37
1532 17545 : t47 = br_c4*t46
1533 17545 : t48 = t39*t29
1534 17545 : t49 = t5*t48
1535 17545 : t50 = t32**2
1536 17545 : t51 = 0.1e1_dp/t50
1537 17545 : t52 = t49*t51
1538 17545 : t56 = t1*t37*pi
1539 17545 : t57 = br_c5*t56
1540 17545 : t58 = t39**2
1541 17545 : t59 = t4*t58
1542 17545 : t61 = 0.1e1_dp/t50/t16
1543 17545 : t62 = t59*t61
1544 : t65 = br_c0 + 0.2e1_dp/0.3e1_dp*t24*t18 + 0.4e1_dp/0.9e1_dp*t28*t34 &
1545 : + 0.8e1_dp/0.27e2_dp*t38*t43 + 0.16e2_dp/0.81e2_dp*t47*t52 + 0.32e2_dp &
1546 17545 : /0.243e3_dp*t57*t62
1547 17545 : t66 = t23*t65
1548 17545 : t67 = br_b1*t2
1549 17545 : t70 = br_b2*t27
1550 17545 : t73 = br_b3*t37
1551 17545 : t76 = br_b4*t46
1552 17545 : t79 = br_b5*t56
1553 : t82 = br_b0 + 0.2e1_dp/0.3e1_dp*t67*t18 + 0.4e1_dp/0.9e1_dp*t70*t34 &
1554 : + 0.8e1_dp/0.27e2_dp*t73*t43 + 0.16e2_dp/0.81e2_dp*t76*t52 + 0.32e2_dp &
1555 17545 : /0.243e3_dp*t79*t62
1556 17545 : t83 = 0.1e1_dp/t82
1557 17545 : t84 = t66*t83
1558 17545 : t85 = 8**(0.1e1_dp/0.3e1_dp)
1559 17545 : t86 = t23**2
1560 17545 : t87 = t86*t23
1561 17545 : t88 = t65**2
1562 17545 : t89 = t88*t65
1563 17545 : t90 = t87*t89
1564 17545 : t91 = t82**2
1565 17545 : t93 = 0.1e1_dp/t91/t82
1566 17545 : t94 = t90*t93
1567 17545 : t95 = EXP(-t84)
1568 17545 : t96 = 0.1e1_dp/0.3141592654e1_dp
1569 17545 : t97 = t95*t96
1570 17545 : t99 = t94*t97*t10
1571 17545 : t100 = t99**(0.1e1_dp/0.3e1_dp)
1572 17545 : t101 = 0.1e1_dp/t100
1573 17545 : t102 = REAL(t85, KIND=dp)*t101
1574 17545 : t103 = t102*R
1575 17545 : t104 = t84*t103
1576 17545 : t106 = EXP(0.2e1_dp*t104)
1577 17545 : t108 = t106*R
1578 17545 : t109 = t108*t23
1579 17545 : t110 = t65*t83
1580 17545 : t111 = t110*t102
1581 17545 : t113 = t106*t23
1582 17545 : t115 = t84 + t104
1583 17545 : t116 = EXP(t115)
1584 : t118 = 0.2e1_dp*t106 - t109*t111 + t113*t110 + 0.2e1_dp + t84 + t104 &
1585 17545 : - 0.4e1_dp*t116
1586 17545 : t119 = rho*t118
1587 17545 : t121 = EXP(-t115)
1588 17545 : t123 = t121*REAL(t85, KIND=dp)*t101
1589 17545 : e_0 = e_0 + (t119*t123/0.8e1_dp)*sx
1590 17545 : IF (grad_deriv >= 1) THEN
1591 9431 : t128 = t5*t17
1592 9431 : t131 = 0.1e1_dp/t4
1593 9431 : t133 = t33*gamma
1594 9431 : t134 = t133*t9
1595 9431 : t137 = 0.10e2_dp/0.9e1_dp*t3*t128 + t3*t131*t134/0.18e2_dp
1596 9431 : t138 = t21**2
1597 9431 : t140 = 0.1e1_dp/(0.1e1_dp + t138)
1598 9431 : t141 = t137*t140
1599 9431 : t143 = t83*REAL(t85, KIND=dp)
1600 9431 : t146 = t141*t65*t143*t101*R
1601 9431 : t153 = t4*t29
1602 9431 : t154 = t153*t33
1603 9431 : t157 = t4*rho
1604 9431 : t159 = t42*gamma
1605 9431 : t160 = t159*t9
1606 9431 : t163 = t39*t42
1607 9431 : t167 = t51*gamma
1608 9431 : t168 = t167*t9
1609 9431 : t171 = t5*t40
1610 9431 : t172 = t171*t51
1611 9431 : t175 = t5*t39
1612 9431 : t177 = t61*gamma
1613 9431 : t178 = t177*t9
1614 9431 : t182 = t4*t39*t30
1615 9431 : t183 = t182*t61
1616 9431 : t186 = t4*t48
1617 9431 : t189 = 0.1e1_dp/t50/t32
1618 9431 : t190 = t189*gamma
1619 9431 : t191 = t190*t9
1620 : t194 = 0.10e2_dp/0.9e1_dp*t24*t128 + t24*t131*t134/0.18e2_dp + 0.40e2_dp &
1621 : /0.27e2_dp*t28*t154 + 0.2e1_dp/0.27e2_dp*t28*t157*t160 + &
1622 : 0.40e2_dp/0.27e2_dp*t38*t163 + 0.2e1_dp/0.27e2_dp*t38*t30*t168 &
1623 : + 0.320e3_dp/0.243e3_dp*t47*t172 + 0.16e2_dp/0.243e3_dp*t47*t175* &
1624 : t178 + 0.800e3_dp/0.729e3_dp*t57*t183 + 0.40e2_dp/0.729e3_dp*t57* &
1625 9431 : t186*t191
1626 9431 : t196 = t23*t194*t83
1627 9431 : t197 = t196*t103
1628 9431 : t199 = 0.1e1_dp/t91
1629 9431 : t200 = t66*t199
1630 : t226 = 0.10e2_dp/0.9e1_dp*t67*t128 + t67*t131*t134/0.18e2_dp + 0.40e2_dp &
1631 : /0.27e2_dp*t70*t154 + 0.2e1_dp/0.27e2_dp*t70*t157*t160 + &
1632 : 0.40e2_dp/0.27e2_dp*t73*t163 + 0.2e1_dp/0.27e2_dp*t73*t30*t168 &
1633 : + 0.320e3_dp/0.243e3_dp*t76*t172 + 0.16e2_dp/0.243e3_dp*t76*t175* &
1634 : t178 + 0.800e3_dp/0.729e3_dp*t79*t183 + 0.40e2_dp/0.729e3_dp*t79* &
1635 9431 : t186*t191
1636 9431 : t229 = t200*t102*R*t226
1637 9431 : t232 = 0.1e1_dp/t100/t99
1638 9431 : t233 = REAL(t85, KIND=dp)*t232
1639 9431 : t234 = t86*t89
1640 9431 : t235 = t93*t95
1641 9431 : t237 = t96*t10
1642 9431 : t242 = t87*t88*t93
1643 9431 : t247 = t91**2
1644 9431 : t249 = t90/t247
1645 9431 : t254 = t141*t110
1646 9431 : t256 = t66*t199*t226
1647 : t264 = -0.3e1_dp*t234*t235*t237*t141 + 0.3e1_dp*t242*t97*t10 &
1648 : *t194 - 0.3e1_dp*t249*t97*t10*t226 + t94*(t254 - t196 + &
1649 9431 : t256)*t95*t237 - t94*t97/t29
1650 9431 : t267 = t84*t233*R*t264
1651 : t270 = (-0.2e1_dp*t146 + 0.2e1_dp*t197 - 0.2e1_dp*t229 - 0.2e1_dp/0.3e1_dp &
1652 9431 : *t267)*t106
1653 9431 : t272 = R*t23
1654 9431 : t277 = t194*t83
1655 9431 : t280 = t108*t66
1656 9431 : t281 = t199*REAL(t85, KIND=dp)
1657 9431 : t292 = t140*t65*t83
1658 9431 : t295 = t65*t199
1659 9431 : t298 = t267/0.3e1_dp
1660 9431 : t299 = -t254 + t196 - t256 - t146 + t197 - t229 - t298
1661 : t302 = 0.2e1_dp*t270 - t270*t272*t111 + t108*t141*t111 - t109 &
1662 : *t277*t102 + t280*t281*t101*t226 + t280*t143*t232* &
1663 : t264/0.3e1_dp + t270*t84 - t106*t137*t292 + t113*t277 - t113 &
1664 : *t295*t226 - t254 + t196 - t256 - t146 + t197 - t229 - t298 &
1665 9431 : - 0.4e1_dp*t299*t116
1666 9431 : t310 = t119*t121
1667 : e_rho = e_rho + (t118*t121*t102/0.8e1_dp + rho*t302*t123/0.8e1_dp - t119 &
1668 9431 : *t299*t123/0.8e1_dp - t310*t233*t264/0.24e2_dp)*sx
1669 9431 : t314 = t3*t5
1670 9431 : t315 = t133*ndrho
1671 9431 : t317 = t292*t103
1672 9431 : t318 = t314*t315*t317
1673 9431 : t324 = t159*ndrho
1674 9431 : t328 = t167*ndrho
1675 9431 : t332 = t177*ndrho
1676 9431 : t336 = t190*ndrho
1677 : t339 = -t24*t5*t315/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t28*t153*t324 &
1678 : - 0.4e1_dp/0.27e2_dp*t38*t39*t328 - 0.32e2_dp/0.243e3_dp*t47*t171 &
1679 9431 : *t332 - 0.80e2_dp/0.729e3_dp*t57*t182*t336
1680 9431 : t341 = t23*t339*t83
1681 9431 : t342 = t341*t103
1682 : t359 = -t67*t5*t315/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t70*t153*t324 &
1683 : - 0.4e1_dp/0.27e2_dp*t73*t39*t328 - 0.32e2_dp/0.243e3_dp*t76*t171 &
1684 9431 : *t332 - 0.80e2_dp/0.729e3_dp*t79*t182*t336
1685 9431 : t362 = t200*t102*R*t359
1686 9431 : t368 = gamma*ndrho
1687 9431 : t369 = t368*t140
1688 9431 : t383 = t368*t292
1689 9431 : t385 = t3*t5*t33*t383/0.9e1_dp
1690 9431 : t387 = t66*t199*t359
1691 : t392 = t234*t93*t97*t131*t3*t33*t369/0.3e1_dp + 0.3e1_dp* &
1692 : t242*t97*t10*t339 - 0.3e1_dp*t249*t97*t10*t359 + t94* &
1693 9431 : (-t385 - t341 + t387)*t95*t237
1694 9431 : t395 = t84*t233*R*t392
1695 : t398 = (0.2e1_dp/0.9e1_dp*t318 + 0.2e1_dp*t342 - 0.2e1_dp*t362 - 0.2e1_dp &
1696 9431 : /0.3e1_dp*t395)*t106
1697 9431 : t402 = t108*br_a1
1698 9431 : t404 = t2*t5*t33
1699 9431 : t409 = t339*t83
1700 9431 : t420 = t106*br_a1
1701 9431 : t427 = t318/0.9e1_dp
1702 9431 : t428 = t395/0.3e1_dp
1703 9431 : t429 = t385 + t341 - t387 + t427 + t342 - t362 - t428
1704 : t432 = 0.2e1_dp*t398 - t398*t272*t111 - t402*t404*t369*t111 &
1705 : /0.9e1_dp - t109*t409*t102 + t280*t281*t101*t359 + t280 &
1706 : *t143*t232*t392/0.3e1_dp + t398*t84 + t420*t404*t383/0.9e1_dp &
1707 : + t113*t409 - t113*t295*t359 + t385 + t341 - t387 + t427 &
1708 9431 : + t342 - t362 - t428 - 0.4e1_dp*t429*t116
1709 : e_ndrho = e_ndrho + (rho*t432*t123/0.8e1_dp - t119*t429*t123/0.8e1_dp - t310 &
1710 9431 : *t233*t392/0.24e2_dp)*sx
1711 9431 : t443 = t6*t33
1712 9431 : t444 = t443*gamma
1713 9431 : t446 = t3*t444*t317
1714 9431 : t450 = t31*t42
1715 9431 : t451 = t450*gamma
1716 9431 : t454 = t40*t51
1717 9431 : t455 = t454*gamma
1718 9431 : t458 = t49*t61
1719 9431 : t459 = t458*gamma
1720 9431 : t462 = t59*t189
1721 9431 : t463 = t462*gamma
1722 : t466 = 0.4e1_dp/0.9e1_dp*t24*t444 + 0.16e2_dp/0.27e2_dp*t28*t451 + &
1723 : 0.16e2_dp/0.27e2_dp*t38*t455 + 0.128e3_dp/0.243e3_dp*t47*t459 + 0.320e3_dp &
1724 9431 : /0.729e3_dp*t57*t463
1725 9431 : t468 = t23*t466*t83
1726 9431 : t469 = t468*t103
1727 : t481 = 0.4e1_dp/0.9e1_dp*t67*t444 + 0.16e2_dp/0.27e2_dp*t70*t451 + &
1728 : 0.16e2_dp/0.27e2_dp*t73*t455 + 0.128e3_dp/0.243e3_dp*t76*t459 + 0.320e3_dp &
1729 9431 : /0.729e3_dp*t79*t463
1730 9431 : t484 = t200*t102*R*t481
1731 9431 : t487 = t234*t235*t96
1732 9431 : t501 = gamma*t140
1733 9431 : t504 = 0.4e1_dp/0.9e1_dp*t3*t443*t501*t110
1734 9431 : t506 = t66*t199*t481
1735 : t511 = -0.4e1_dp/0.3e1_dp*t487*t314*t133*t140 + 0.3e1_dp*t242* &
1736 : t97*t10*t466 - 0.3e1_dp*t249*t97*t10*t481 + t94*(t504 - &
1737 9431 : t468 + t506)*t95*t237
1738 9431 : t514 = t84*t233*R*t511
1739 : t517 = (-0.8e1_dp/0.9e1_dp*t446 + 0.2e1_dp*t469 - 0.2e1_dp*t484 - 0.2e1_dp &
1740 9431 : /0.3e1_dp*t514)*t106
1741 9431 : t521 = t2*t6
1742 9431 : t525 = t143*t101
1743 9431 : t529 = t466*t83
1744 9431 : t540 = t420*t521
1745 9431 : t547 = 0.4e1_dp/0.9e1_dp*t446
1746 9431 : t548 = t514/0.3e1_dp
1747 9431 : t549 = -t504 + t468 - t506 - t547 + t469 - t484 - t548
1748 : t552 = 0.2e1_dp*t517 - t517*t272*t111 + 0.4e1_dp/0.9e1_dp*t402*t521 &
1749 : *t33*t501*t65*t525 - t109*t529*t102 + t280*t281* &
1750 : t101*t481 + t280*t143*t232*t511/0.3e1_dp + t517*t84 - 0.4e1_dp &
1751 : /0.9e1_dp*t540*t133*t292 + t113*t529 - t113*t295*t481 &
1752 : - t504 + t468 - t506 - t547 + t469 - t484 - t548 - 0.4e1_dp*t549 &
1753 9431 : *t116
1754 : e_tau = e_tau + (rho*t552*t123/0.8e1_dp - t119*t549*t123/0.8e1_dp - t310 &
1755 9431 : *t233*t511/0.24e2_dp)*sx
1756 9431 : t566 = t3*t443*t140*t110*t103
1757 : t578 = -t24*t443/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t28*t450 - 0.4e1_dp/ &
1758 : 0.27e2_dp*t38*t454 - 0.32e2_dp/0.243e3_dp*t47*t458 - 0.80e2_dp/0.729e3_dp &
1759 9431 : *t57*t462
1760 9431 : t580 = t23*t578*t83
1761 9431 : t581 = t580*t103
1762 : t593 = -t67*t443/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t70*t450 - 0.4e1_dp/ &
1763 : 0.27e2_dp*t73*t454 - 0.32e2_dp/0.243e3_dp*t76*t458 - 0.80e2_dp/0.729e3_dp &
1764 9431 : *t79*t462
1765 9431 : t596 = t200*t102*R*t593
1766 9431 : t612 = t3*t6
1767 9431 : t613 = t33*t140
1768 9431 : t614 = t613*t110
1769 9431 : t616 = t612*t614/0.9e1_dp
1770 9431 : t618 = t66*t199*t593
1771 : t623 = t487*t5*br_a1*t2*t33*t140/0.3e1_dp + 0.3e1_dp*t242* &
1772 : t97*t10*t578 - 0.3e1_dp*t249*t97*t10*t593 + t94*(-t616 &
1773 9431 : - t580 + t618)*t95*t237
1774 9431 : t626 = t84*t233*R*t623
1775 : t629 = (0.2e1_dp/0.9e1_dp*t566 + 0.2e1_dp*t581 - 0.2e1_dp*t596 - 0.2e1_dp &
1776 9431 : /0.3e1_dp*t626)*t106
1777 9431 : t638 = t578*t83
1778 9431 : t654 = t566/0.9e1_dp
1779 9431 : t655 = t626/0.3e1_dp
1780 9431 : t656 = t616 + t580 - t618 + t654 + t581 - t596 - t655
1781 : t659 = 0.2e1_dp*t629 - t629*t272*t111 - t108*t612*t613*t65 &
1782 : *t525/0.9e1_dp - t109*t638*t102 + t280*t281*t101*t593 + &
1783 : t280*t143*t232*t623/0.3e1_dp + t629*t84 + t540*t614/0.9e1_dp &
1784 : + t113*t638 - t113*t295*t593 + t616 + t580 - t618 + t654 &
1785 9431 : + t581 - t596 - t655 - 0.4e1_dp*t656*t116
1786 : e_laplace_rho = e_laplace_rho + (rho*t659*t123/0.8e1_dp - t119*t656*t123/0.8e1_dp - t310 &
1787 9431 : *t233*t623/0.24e2_dp)*sx
1788 : END IF
1789 :
1790 17545 : END SUBROUTINE x_br_lsd_y_lte_0_cutoff_R_lte_b
1791 :
1792 : ! **************************************************************************************************
1793 : !> \brief Truncated long range part for y > 0
1794 : !> \param rho ...
1795 : !> \param ndrho ...
1796 : !> \param tau ...
1797 : !> \param laplace_rho ...
1798 : !> \param e_0 ...
1799 : !> \param e_rho ...
1800 : !> \param e_ndrho ...
1801 : !> \param e_tau ...
1802 : !> \param e_laplace_rho ...
1803 : !> \param sx ...
1804 : !> \param R ...
1805 : !> \param gamma ...
1806 : !> \param grad_deriv ...
1807 : !> \par History
1808 : !> 12.2008 created [mguidon]
1809 : !> \author mguidon
1810 : ! **************************************************************************************************
1811 16977565 : SUBROUTINE x_br_lsd_y_gt_0_cutoff(rho, ndrho, tau, laplace_rho, e_0, &
1812 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1813 : sx, R, gamma, grad_deriv)
1814 : REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho
1815 : REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho
1816 : REAL(dp), INTENT(IN) :: sx, R, gamma
1817 : INTEGER, INTENT(IN) :: grad_deriv
1818 :
1819 : REAL(KIND=dp) :: brval, t1, t10, t103, t104, t11, t111, t116, t14, t15, t2, t21, t25, t26, &
1820 : t28, t3, t31, t33, t37, t4, t44, t45, t46, t5, t50, t56, t58, t6, t62, t65, t69, t71, &
1821 : t75, t77, t8, t81, t84, t85, t9
1822 :
1823 16977565 : t1 = 8**(0.1e1_dp/0.3e1_dp)
1824 16977565 : t2 = t1**2
1825 16977565 : t3 = 0.1e1_dp/br_BB
1826 16977565 : t4 = pi**(0.1e1_dp/0.3e1_dp)
1827 16977565 : t5 = t4**2
1828 16977565 : t6 = 0.1e1_dp/t5
1829 16977565 : t8 = rho**(0.1e1_dp/0.3e1_dp)
1830 16977565 : t9 = t8**2
1831 16977565 : t10 = t9*rho
1832 16977565 : t11 = 0.1e1_dp/t10
1833 16977565 : t14 = ndrho**2
1834 16977565 : t15 = 0.1e1_dp/rho
1835 16977565 : t21 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t14*t15/0.4e1_dp)/0.3e1_dp
1836 16977565 : t25 = br_BB**2
1837 16977565 : t26 = t4*pi
1838 16977565 : t28 = rho**2
1839 16977565 : t31 = t21**2
1840 16977565 : t33 = t8*t28*rho/t31
1841 16977565 : t37 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t26*t33)
1842 : t44 = LOG(0.1500000000000000e1_dp*t3*t6*t11*t21 + 0.3e1_dp/0.2e1_dp &
1843 16977565 : *t37*t3*t6*t11*t21)
1844 16977565 : t45 = t44 + 0.2e1_dp
1845 16977565 : t46 = t45**2
1846 16977565 : t50 = t10/t21
1847 16977565 : t56 = pi**2
1848 16977565 : t58 = t28**2
1849 16977565 : t62 = t58*rho/t31/t21
1850 16977565 : t65 = t5*t56
1851 16977565 : t69 = t31**2
1852 16977565 : t71 = t9*t58*t28/t69
1853 16977565 : t75 = t4*t56*pi
1854 16977565 : t77 = t58**2
1855 16977565 : t81 = t8*t77/t69/t21
1856 : t84 = br_d0 + 0.2e1_dp/0.3e1_dp*br_d1*t5*t50 + 0.4e1_dp/0.9e1_dp*br_d2 &
1857 : *t26*t33 + 0.8e1_dp/0.27e2_dp*br_d3*t56*t62 + 0.16e2_dp/0.81e2_dp &
1858 16977565 : *br_d4*t65*t71 + 0.32e2_dp/0.243e3_dp*br_d5*t75*t81
1859 16977565 : t85 = t84**2
1860 : t103 = br_e0 + 0.2e1_dp/0.3e1_dp*br_e1*t5*t50 + 0.4e1_dp/0.9e1_dp*br_e2 &
1861 : *t26*t33 + 0.8e1_dp/0.27e2_dp*br_e3*t56*t62 + 0.16e2_dp/0.81e2_dp &
1862 16977565 : *br_e4*t65*t71 + 0.32e2_dp/0.243e3_dp*br_e5*t75*t81
1863 16977565 : t104 = t103**2
1864 16977565 : t111 = EXP(-t45*t84/t103)
1865 : t116 = (t46*t45*t85*t84/t104/t103*t111/0.3141592654e1_dp &
1866 16977565 : *t15)**(0.1e1_dp/0.3e1_dp)
1867 16977565 : brval = REAL(t2, KIND=dp)*t116/0.8e1_dp
1868 :
1869 16977565 : IF (R > brval) THEN
1870 : CALL x_br_lsd_y_gt_0_cutoff_R_gt_b(rho, ndrho, tau, laplace_rho, e_0, &
1871 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1872 804729 : sx, R, gamma, grad_deriv)
1873 : ELSE
1874 : CALL x_br_lsd_y_gt_0_cutoff_R_lte_b(rho, ndrho, tau, laplace_rho, e_0, &
1875 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1876 16172836 : sx, R, gamma, grad_deriv)
1877 : END IF
1878 :
1879 16977565 : END SUBROUTINE x_br_lsd_y_gt_0_cutoff
1880 :
1881 : ! **************************************************************************************************
1882 : !> \brief Truncated long range part for y > 0 and R > b
1883 : !> \param rho ...
1884 : !> \param ndrho ...
1885 : !> \param tau ...
1886 : !> \param laplace_rho ...
1887 : !> \param e_0 ...
1888 : !> \param e_rho ...
1889 : !> \param e_ndrho ...
1890 : !> \param e_tau ...
1891 : !> \param e_laplace_rho ...
1892 : !> \param sx ...
1893 : !> \param R ...
1894 : !> \param gamma ...
1895 : !> \param grad_deriv ...
1896 : !> \par History
1897 : !> 12.2008 created [mguidon]
1898 : !> \author mguidon
1899 : ! **************************************************************************************************
1900 804729 : SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_gt_b(rho, ndrho, tau, laplace_rho, e_0, &
1901 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
1902 : sx, R, gamma, grad_deriv)
1903 :
1904 : REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho
1905 : REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho
1906 : REAL(dp), INTENT(IN) :: sx, R, gamma
1907 : INTEGER, INTENT(IN) :: grad_deriv
1908 :
1909 : REAL(KIND=dp) :: t1, t100, t101, t102, t103, t104, t105, t106, t108, t109, t110, t111, t112, &
1910 : t114, t115, t116, t117, t118, t119, t12, t121, t123, t124, t125, t129, t13, t132, t134, &
1911 : t135, t138, t139, t145, t152, t155, t158, t159, t162, t164, t165, t176, t179, t180, t181, &
1912 : t182, t183, t186, t188, t189, t19, t197, t2, t20, t201, t202, t205, t206, t209, t211, &
1913 : t212, t216, t217, t220, t223, t224, t225, t228, t23, t230, t231, t24, t25, t257, t259, &
1914 : t26, t261, t262, t263, t265, t266, t269, t27, t272, t273, t278, t28, t283, t285, t29, &
1915 : t297, t299, t3, t30, t301, t303, t306, t307, t308, t31, t310, t312
1916 : REAL(KIND=dp) :: t313, t321, t326, t329, t339, t343, t35, t351, t354, t355, t36, t363, t364, &
1917 : t365, t367, t37, t371, t375, t379, t383, t386, t388, t4, t404, t406, t408, t409, t41, &
1918 : t411, t412, t42, t428, t43, t430, t432, t434, t437, t439, t44, t441, t45, t453, t456, &
1919 : t46, t469, t479, t480, t485, t486, t487, t49, t490, t491, t494, t495, t498, t499, t5, &
1920 : t502, t503, t506, t508, t519, t52, t521, t523, t524, t526, t527, t53, t54, t543, t545, &
1921 : t547, t549, t55, t552, t554, t556, t568, t57, t571, t58, t584, t599, t6, t600, t601, t61, &
1922 : t612, t614, t62, t625, t627, t629, t63, t630, t632, t633, t64, t649
1923 : REAL(KIND=dp) :: t65, t651, t653, t655, t658, t66, t660, t662, t67, t674, t677, t690, t7, &
1924 : t71, t72, t73, t74, t76, t77, t8, t80, t81, t82, t85, t88, t9, t91, t94, t97, t98, t99
1925 :
1926 804729 : t1 = 0.1e1_dp/br_BB
1927 804729 : t2 = pi**(0.1e1_dp/0.3e1_dp)
1928 804729 : t3 = t2**2
1929 804729 : t4 = 0.1e1_dp/t3
1930 804729 : t5 = t1*t4
1931 804729 : t6 = rho**(0.1e1_dp/0.3e1_dp)
1932 804729 : t7 = t6**2
1933 804729 : t8 = t7*rho
1934 804729 : t9 = 0.1e1_dp/t8
1935 804729 : t12 = ndrho**2
1936 804729 : t13 = 0.1e1_dp/rho
1937 804729 : t19 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t12*t13/0.4e1_dp)/0.3e1_dp
1938 804729 : t20 = t9*t19
1939 804729 : t23 = br_BB**2
1940 804729 : t24 = t2*pi
1941 804729 : t25 = t23*t24
1942 804729 : t26 = rho**2
1943 804729 : t27 = t26*rho
1944 804729 : t28 = t6*t27
1945 804729 : t29 = t19**2
1946 804729 : t30 = 0.1e1_dp/t29
1947 804729 : t31 = t28*t30
1948 804729 : t35 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t31)
1949 804729 : t36 = t35*t1
1950 804729 : t37 = t4*t9
1951 : t41 = 0.1500000000000000e1_dp*t5*t20 + 0.3e1_dp/0.2e1_dp*t36*t37* &
1952 804729 : t19
1953 804729 : t42 = LOG(t41)
1954 804729 : t43 = t42 + 0.2e1_dp
1955 804729 : t44 = br_d1*t3
1956 804729 : t45 = 0.1e1_dp/t19
1957 804729 : t46 = t8*t45
1958 804729 : t49 = br_d2*t24
1959 804729 : t52 = pi**2
1960 804729 : t53 = br_d3*t52
1961 804729 : t54 = t26**2
1962 804729 : t55 = t54*rho
1963 804729 : t57 = 0.1e1_dp/t29/t19
1964 804729 : t58 = t55*t57
1965 804729 : t61 = t3*t52
1966 804729 : t62 = br_d4*t61
1967 804729 : t63 = t54*t26
1968 804729 : t64 = t7*t63
1969 804729 : t65 = t29**2
1970 804729 : t66 = 0.1e1_dp/t65
1971 804729 : t67 = t64*t66
1972 804729 : t71 = t2*t52*pi
1973 804729 : t72 = br_d5*t71
1974 804729 : t73 = t54**2
1975 804729 : t74 = t6*t73
1976 804729 : t76 = 0.1e1_dp/t65/t19
1977 804729 : t77 = t74*t76
1978 : t80 = br_d0 + 0.2e1_dp/0.3e1_dp*t44*t46 + 0.4e1_dp/0.9e1_dp*t49*t31 &
1979 : + 0.8e1_dp/0.27e2_dp*t53*t58 + 0.16e2_dp/0.81e2_dp*t62*t67 + 0.32e2_dp &
1980 804729 : /0.243e3_dp*t72*t77
1981 804729 : t81 = t43*t80
1982 804729 : t82 = br_e1*t3
1983 804729 : t85 = br_e2*t24
1984 804729 : t88 = br_e3*t52
1985 804729 : t91 = br_e4*t61
1986 804729 : t94 = br_e5*t71
1987 : t97 = br_e0 + 0.2e1_dp/0.3e1_dp*t82*t46 + 0.4e1_dp/0.9e1_dp*t85*t31 &
1988 : + 0.8e1_dp/0.27e2_dp*t88*t58 + 0.16e2_dp/0.81e2_dp*t91*t67 + 0.32e2_dp &
1989 804729 : /0.243e3_dp*t94*t77
1990 804729 : t98 = 0.1e1_dp/t97
1991 804729 : t99 = t81*t98
1992 804729 : t100 = 8**(0.1e1_dp/0.3e1_dp)
1993 804729 : t101 = t43**2
1994 804729 : t102 = t101*t43
1995 804729 : t103 = t80**2
1996 804729 : t104 = t103*t80
1997 804729 : t105 = t102*t104
1998 804729 : t106 = t97**2
1999 804729 : t108 = 0.1e1_dp/t106/t97
2000 804729 : t109 = t105*t108
2001 804729 : t110 = EXP(-t99)
2002 804729 : t111 = 0.1e1_dp/0.3141592654e1_dp
2003 804729 : t112 = t110*t111
2004 804729 : t114 = t109*t112*t13
2005 804729 : t115 = t114**(0.1e1_dp/0.3e1_dp)
2006 804729 : t116 = 0.1e1_dp/t115
2007 804729 : t117 = REAL(t100, KIND=dp)*t116
2008 804729 : t118 = t117*R
2009 804729 : t119 = t99*t118
2010 804729 : t121 = EXP(t99 - t119)
2011 804729 : t123 = t121*t43
2012 804729 : t124 = t80*t98
2013 804729 : t125 = t123*t124
2014 804729 : t129 = t98*REAL(t100, KIND=dp)*t116*R
2015 804729 : t132 = EXP(-t99 - t119)
2016 804729 : t134 = t132*t43
2017 804729 : t135 = t134*t124
2018 : t138 = -0.2e1_dp*t121 + t125 - t123*t80*t129 + 0.2e1_dp*t132 + t135 &
2019 804729 : + t134*t80*t129
2020 804729 : t139 = rho*t138
2021 804729 : e_0 = e_0 + (t139*t117/0.8e1_dp)*sx
2022 804729 : IF (grad_deriv >= 1) THEN
2023 415000 : t145 = 0.1e1_dp/t7/t26
2024 415000 : t152 = 0.1e1_dp/t7/t27*gamma*t12
2025 415000 : t155 = 0.1e1_dp/t35
2026 415000 : t158 = t6*t26
2027 415000 : t159 = t158*t30
2028 415000 : t162 = t6*rho
2029 415000 : t164 = t57*gamma
2030 415000 : t165 = t164*t12
2031 415000 : t176 = t36*t4
2032 : t179 = -0.2500000000e1_dp*t5*t145*t19 - 0.1250000000e0_dp*t5*t152 &
2033 : + 0.3e1_dp/0.4e1_dp*t155*t1*t4*t20*(0.40e2_dp/0.27e2_dp*t25 &
2034 : *t159 + 0.2e1_dp/0.27e2_dp*t25*t162*t165) - 0.5e1_dp/0.2e1_dp*t36 &
2035 415000 : *t4*t145*t19 - t176*t152/0.8e1_dp
2036 415000 : t180 = 0.1e1_dp/t41
2037 415000 : t181 = t179*t180
2038 415000 : t182 = t181*t124
2039 415000 : t183 = t7*t45
2040 415000 : t186 = 0.1e1_dp/t6
2041 415000 : t188 = t30*gamma
2042 415000 : t189 = t188*t12
2043 415000 : t197 = t54*t57
2044 415000 : t201 = t66*gamma
2045 415000 : t202 = t201*t12
2046 415000 : t205 = t7*t55
2047 415000 : t206 = t205*t66
2048 415000 : t209 = t7*t54
2049 415000 : t211 = t76*gamma
2050 415000 : t212 = t211*t12
2051 415000 : t216 = t6*t54*t27
2052 415000 : t217 = t216*t76
2053 415000 : t220 = t6*t63
2054 415000 : t223 = 0.1e1_dp/t65/t29
2055 415000 : t224 = t223*gamma
2056 415000 : t225 = t224*t12
2057 : t228 = 0.10e2_dp/0.9e1_dp*t44*t183 + t44*t186*t189/0.18e2_dp + 0.40e2_dp &
2058 : /0.27e2_dp*t49*t159 + 0.2e1_dp/0.27e2_dp*t49*t162*t165 + &
2059 : 0.40e2_dp/0.27e2_dp*t53*t197 + 0.2e1_dp/0.27e2_dp*t53*t27*t202 &
2060 : + 0.320e3_dp/0.243e3_dp*t62*t206 + 0.16e2_dp/0.243e3_dp*t62*t209* &
2061 : t212 + 0.800e3_dp/0.729e3_dp*t72*t217 + 0.40e2_dp/0.729e3_dp*t72* &
2062 415000 : t220*t225
2063 415000 : t230 = t43*t228*t98
2064 415000 : t231 = 0.1e1_dp/t106
2065 : t257 = 0.10e2_dp/0.9e1_dp*t82*t183 + t82*t186*t189/0.18e2_dp + 0.40e2_dp &
2066 : /0.27e2_dp*t85*t159 + 0.2e1_dp/0.27e2_dp*t85*t162*t165 + &
2067 : 0.40e2_dp/0.27e2_dp*t88*t197 + 0.2e1_dp/0.27e2_dp*t88*t27*t202 &
2068 : + 0.320e3_dp/0.243e3_dp*t91*t206 + 0.16e2_dp/0.243e3_dp*t91*t209* &
2069 : t212 + 0.800e3_dp/0.729e3_dp*t94*t217 + 0.40e2_dp/0.729e3_dp*t94* &
2070 415000 : t220*t225
2071 415000 : t259 = t81*t231*t257
2072 415000 : t261 = t181*t80*t129
2073 415000 : t262 = t230*t118
2074 415000 : t263 = t81*t231
2075 415000 : t265 = t117*R*t257
2076 415000 : t266 = t263*t265
2077 415000 : t269 = REAL(t100, KIND=dp)/t115/t114
2078 415000 : t272 = t101*t104*t108*t110
2079 415000 : t273 = t111*t13
2080 415000 : t278 = t102*t103*t108
2081 415000 : t283 = t106**2
2082 415000 : t285 = t105/t283
2083 : t297 = 0.3e1_dp*t272*t273*t181 + 0.3e1_dp*t278*t112*t13*t228 &
2084 : - 0.3e1_dp*t285*t112*t13*t257 + t109*(-t182 - t230 + t259) &
2085 415000 : *t110*t273 - t109*t112/t26
2086 415000 : t299 = t269*R*t297
2087 415000 : t301 = t99*t299/0.3e1_dp
2088 415000 : t303 = (t182 + t230 - t259 - t261 - t262 + t266 + t301)*t121
2089 415000 : t306 = t121*t179
2090 415000 : t307 = t180*t80
2091 415000 : t308 = t307*t98
2092 415000 : t310 = t228*t98
2093 415000 : t312 = t80*t231
2094 415000 : t313 = t312*t257
2095 415000 : t321 = t123*t312
2096 415000 : t326 = (-t182 - t230 + t259 - t261 - t262 + t266 + t301)*t132
2097 415000 : t329 = t132*t179
2098 415000 : t339 = t134*t312
2099 : t343 = -0.2e1_dp*t303 + t303*t99 + t306*t308 + t123*t310 - t123 &
2100 : *t313 - t303*t81*t129 - t306*t307*t129 - t123*t228* &
2101 : t129 + t321*t265 + t125*t299/0.3e1_dp + 0.2e1_dp*t326 + t326*t99 &
2102 : + t329*t308 + t134*t310 - t134*t313 + t326*t81*t129 + &
2103 : t329*t307*t129 + t134*t228*t129 - t339*t265 - t135*t299 &
2104 415000 : /0.3e1_dp
2105 : e_rho = e_rho + (t138*REAL(t100, KIND=dp)*t116/0.8e1_dp + rho*t343*t117/0.8e1_dp &
2106 415000 : - t139*t269*t297/0.24e2_dp)*sx
2107 415000 : t351 = t145*gamma*ndrho
2108 415000 : t354 = t155*br_BB
2109 415000 : t355 = t354*t3
2110 : t363 = 0.2500000000000000e0_dp*t5*t351 - t355*t7*t30*gamma*ndrho &
2111 415000 : /0.9e1_dp + t176*t351/0.4e1_dp
2112 415000 : t364 = t363*t180
2113 415000 : t365 = t364*t124
2114 415000 : t367 = t188*ndrho
2115 415000 : t371 = t164*ndrho
2116 415000 : t375 = t201*ndrho
2117 415000 : t379 = t211*ndrho
2118 415000 : t383 = t224*ndrho
2119 : t386 = -t44*t7*t367/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t49*t158*t371 &
2120 : - 0.4e1_dp/0.27e2_dp*t53*t54*t375 - 0.32e2_dp/0.243e3_dp*t62*t205 &
2121 415000 : *t379 - 0.80e2_dp/0.729e3_dp*t72*t216*t383
2122 415000 : t388 = t43*t386*t98
2123 : t404 = -t82*t7*t367/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t85*t158*t371 &
2124 : - 0.4e1_dp/0.27e2_dp*t88*t54*t375 - 0.32e2_dp/0.243e3_dp*t91*t205 &
2125 415000 : *t379 - 0.80e2_dp/0.729e3_dp*t94*t216*t383
2126 415000 : t406 = t81*t231*t404
2127 415000 : t408 = t364*t80*t129
2128 415000 : t409 = t388*t118
2129 415000 : t411 = t117*R*t404
2130 415000 : t412 = t263*t411
2131 : t428 = 0.3e1_dp*t272*t273*t364 + 0.3e1_dp*t278*t112*t13*t386 &
2132 : - 0.3e1_dp*t285*t112*t13*t404 + t109*(-t365 - t388 + t406) &
2133 415000 : *t110*t273
2134 415000 : t430 = t269*R*t428
2135 415000 : t432 = t99*t430/0.3e1_dp
2136 415000 : t434 = (t365 + t388 - t406 - t408 - t409 + t412 + t432)*t121
2137 415000 : t437 = t121*t363
2138 415000 : t439 = t386*t98
2139 415000 : t441 = t312*t404
2140 415000 : t453 = (-t365 - t388 + t406 - t408 - t409 + t412 + t432)*t132
2141 415000 : t456 = t132*t363
2142 : t469 = -0.2e1_dp*t434 + t434*t99 + t437*t308 + t123*t439 - t123 &
2143 : *t441 - t434*t81*t129 - t437*t307*t129 - t123*t386* &
2144 : t129 + t321*t411 + t125*t430/0.3e1_dp + 0.2e1_dp*t453 + t453*t99 &
2145 : + t456*t308 + t134*t439 - t134*t441 + t453*t81*t129 + &
2146 : t456*t307*t129 + t134*t386*t129 - t339*t411 - t135*t430 &
2147 415000 : /0.3e1_dp
2148 415000 : e_ndrho = e_ndrho + (rho*t469*t117/0.8e1_dp - t139*t269*t428/0.24e2_dp)*sx
2149 415000 : t479 = t8*t30
2150 415000 : t480 = t479*gamma
2151 : t485 = -0.1000000000e1_dp*t5*t9*gamma + 0.4e1_dp/0.9e1_dp*t355*t480 &
2152 415000 : - t36*t37*gamma
2153 415000 : t486 = t485*t180
2154 415000 : t487 = t486*t124
2155 415000 : t490 = t28*t57
2156 415000 : t491 = t490*gamma
2157 415000 : t494 = t55*t66
2158 415000 : t495 = t494*gamma
2159 415000 : t498 = t64*t76
2160 415000 : t499 = t498*gamma
2161 415000 : t502 = t74*t223
2162 415000 : t503 = t502*gamma
2163 : t506 = 0.4e1_dp/0.9e1_dp*t44*t480 + 0.16e2_dp/0.27e2_dp*t49*t491 + &
2164 : 0.16e2_dp/0.27e2_dp*t53*t495 + 0.128e3_dp/0.243e3_dp*t62*t499 + 0.320e3_dp &
2165 415000 : /0.729e3_dp*t72*t503
2166 415000 : t508 = t43*t506*t98
2167 : t519 = 0.4e1_dp/0.9e1_dp*t82*t480 + 0.16e2_dp/0.27e2_dp*t85*t491 + &
2168 : 0.16e2_dp/0.27e2_dp*t88*t495 + 0.128e3_dp/0.243e3_dp*t91*t499 + 0.320e3_dp &
2169 415000 : /0.729e3_dp*t94*t503
2170 415000 : t521 = t81*t231*t519
2171 415000 : t523 = t486*t80*t129
2172 415000 : t524 = t508*t118
2173 415000 : t526 = t117*R*t519
2174 415000 : t527 = t263*t526
2175 : t543 = 0.3e1_dp*t272*t273*t486 + 0.3e1_dp*t278*t112*t13*t506 &
2176 : - 0.3e1_dp*t285*t112*t13*t519 + t109*(-t487 - t508 + t521) &
2177 415000 : *t110*t273
2178 415000 : t545 = t269*R*t543
2179 415000 : t547 = t99*t545/0.3e1_dp
2180 415000 : t549 = (t487 + t508 - t521 - t523 - t524 + t527 + t547)*t121
2181 415000 : t552 = t121*t485
2182 415000 : t554 = t506*t98
2183 415000 : t556 = t312*t519
2184 415000 : t568 = (-t487 - t508 + t521 - t523 - t524 + t527 + t547)*t132
2185 415000 : t571 = t132*t485
2186 : t584 = -0.2e1_dp*t549 + t549*t99 + t552*t308 + t123*t554 - t123 &
2187 : *t556 - t549*t81*t129 - t552*t307*t129 - t123*t506* &
2188 : t129 + t321*t526 + t125*t545/0.3e1_dp + 0.2e1_dp*t568 + t568*t99 &
2189 : + t571*t308 + t134*t554 - t134*t556 + t568*t81*t129 + &
2190 : t571*t307*t129 + t134*t506*t129 - t339*t526 - t135*t545 &
2191 415000 : /0.3e1_dp
2192 415000 : e_tau = e_tau + (rho*t584*t117/0.8e1_dp - t139*t269*t543/0.24e2_dp)*sx
2193 : t599 = 0.2500000000000000e0_dp*t5*t9 - t354*t3*t8*t30/0.9e1_dp &
2194 415000 : + t36*t37/0.4e1_dp
2195 415000 : t600 = t599*t180
2196 415000 : t601 = t600*t124
2197 : t612 = -t44*t479/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t49*t490 - 0.4e1_dp/ &
2198 : 0.27e2_dp*t53*t494 - 0.32e2_dp/0.243e3_dp*t62*t498 - 0.80e2_dp/0.729e3_dp &
2199 415000 : *t72*t502
2200 415000 : t614 = t43*t612*t98
2201 : t625 = -t82*t479/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t85*t490 - 0.4e1_dp/ &
2202 : 0.27e2_dp*t88*t494 - 0.32e2_dp/0.243e3_dp*t91*t498 - 0.80e2_dp/0.729e3_dp &
2203 415000 : *t94*t502
2204 415000 : t627 = t81*t231*t625
2205 415000 : t629 = t600*t80*t129
2206 415000 : t630 = t614*t118
2207 415000 : t632 = t117*R*t625
2208 415000 : t633 = t263*t632
2209 : t649 = 0.3e1_dp*t272*t273*t600 + 0.3e1_dp*t278*t112*t13*t612 &
2210 : - 0.3e1_dp*t285*t112*t13*t625 + t109*(-t601 - t614 + t627) &
2211 415000 : *t110*t273
2212 415000 : t651 = t269*R*t649
2213 415000 : t653 = t99*t651/0.3e1_dp
2214 415000 : t655 = (t601 + t614 - t627 - t629 - t630 + t633 + t653)*t121
2215 415000 : t658 = t121*t599
2216 415000 : t660 = t612*t98
2217 415000 : t662 = t312*t625
2218 415000 : t674 = (-t601 - t614 + t627 - t629 - t630 + t633 + t653)*t132
2219 415000 : t677 = t132*t599
2220 : t690 = -0.2e1_dp*t655 + t655*t99 + t658*t308 + t123*t660 - t123 &
2221 : *t662 - t655*t81*t129 - t658*t307*t129 - t123*t612* &
2222 : t129 + t321*t632 + t125*t651/0.3e1_dp + 0.2e1_dp*t674 + t674*t99 &
2223 : + t677*t308 + t134*t660 - t134*t662 + t674*t81*t129 + &
2224 : t677*t307*t129 + t134*t612*t129 - t339*t632 - t135*t651 &
2225 415000 : /0.3e1_dp
2226 415000 : e_laplace_rho = e_laplace_rho + (rho*t690*t117/0.8e1_dp - t139*t269*t649/0.24e2_dp)*sx
2227 : END IF
2228 :
2229 804729 : END SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_gt_b
2230 :
2231 : ! **************************************************************************************************
2232 : !> \brief Truncated long range part for y > 0 and R <= b
2233 : !> \param rho ...
2234 : !> \param ndrho ...
2235 : !> \param tau ...
2236 : !> \param laplace_rho ...
2237 : !> \param e_0 ...
2238 : !> \param e_rho ...
2239 : !> \param e_ndrho ...
2240 : !> \param e_tau ...
2241 : !> \param e_laplace_rho ...
2242 : !> \param sx ...
2243 : !> \param R ...
2244 : !> \param gamma ...
2245 : !> \param grad_deriv ...
2246 : !> \par History
2247 : !> 12.2008 created [mguidon]
2248 : !> \author mguidon
2249 : ! **************************************************************************************************
2250 16172836 : SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_lte_b(rho, ndrho, tau, laplace_rho, e_0, &
2251 : e_rho, e_ndrho, e_tau, e_laplace_rho, &
2252 : sx, R, gamma, grad_deriv)
2253 : REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho
2254 : REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho
2255 : REAL(dp), INTENT(IN) :: sx, R, gamma
2256 : INTEGER, INTENT(IN) :: grad_deriv
2257 :
2258 : REAL(KIND=dp) :: t1, t100, t101, t102, t103, t104, t105, t106, t108, t109, t110, t111, t112, &
2259 : t114, t115, t116, t117, t118, t119, t12, t121, t123, t124, t125, t126, t128, t13, t130, &
2260 : t131, t133, t134, t136, t138, t144, t151, t154, t157, t158, t161, t163, t164, t175, t178, &
2261 : t179, t180, t182, t184, t185, t187, t19, t190, t192, t193, t2, t20, t201, t205, t206, &
2262 : t209, t210, t213, t215, t216, t220, t221, t224, t227, t228, t229, t23, t232, t234, t235, &
2263 : t237, t238, t24, t25, t26, t264, t267, t27, t270, t271, t274, t275, t28, t280, t285, &
2264 : t287, t29, t292, t294, t3, t30, t302, t305, t308, t31, t310, t315
2265 : REAL(KIND=dp) :: t318, t319, t330, t333, t336, t337, t340, t348, t35, t353, t356, t357, t36, &
2266 : t365, t366, t368, t37, t371, t375, t379, t383, t387, t390, t392, t393, t4, t41, t410, &
2267 : t413, t42, t426, t428, t43, t433, t436, t439, t44, t445, t45, t46, t461, t462, t465, &
2268 : t479, t480, t485, t486, t488, t49, t492, t493, t496, t497, t5, t500, t501, t504, t505, &
2269 : t508, t510, t511, t52, t523, t526, t53, t539, t54, t541, t546, t549, t55, t552, t558, &
2270 : t57, t574, t575, t578, t58, t597, t598, t6, t600, t61, t612, t614, t615, t62, t627, t63, &
2271 : t630, t64, t643, t645, t65, t650, t653, t656, t66, t662, t67, t678
2272 : REAL(KIND=dp) :: t679, t682, t7, t71, t72, t73, t74, t76, t77, t8, t80, t81, t82, t85, t88, &
2273 : t9, t91, t94, t97, t98, t99
2274 :
2275 16172836 : t1 = 0.1e1_dp/br_BB
2276 16172836 : t2 = pi**(0.1e1_dp/0.3e1_dp)
2277 16172836 : t3 = t2**2
2278 16172836 : t4 = 0.1e1_dp/t3
2279 16172836 : t5 = t1*t4
2280 16172836 : t6 = rho**(0.1e1_dp/0.3e1_dp)
2281 16172836 : t7 = t6**2
2282 16172836 : t8 = t7*rho
2283 16172836 : t9 = 0.1e1_dp/t8
2284 16172836 : t12 = ndrho**2
2285 16172836 : t13 = 0.1e1_dp/rho
2286 16172836 : t19 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t12*t13/0.4e1_dp)/0.3e1_dp
2287 16172836 : t20 = t9*t19
2288 16172836 : t23 = br_BB**2
2289 16172836 : t24 = t2*pi
2290 16172836 : t25 = t23*t24
2291 16172836 : t26 = rho**2
2292 16172836 : t27 = t26*rho
2293 16172836 : t28 = t6*t27
2294 16172836 : t29 = t19**2
2295 16172836 : t30 = 0.1e1_dp/t29
2296 16172836 : t31 = t28*t30
2297 16172836 : t35 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t31)
2298 16172836 : t36 = t35*t1
2299 16172836 : t37 = t4*t9
2300 : t41 = 0.1500000000000000e1_dp*t5*t20 + 0.3e1_dp/0.2e1_dp*t36*t37* &
2301 16172836 : t19
2302 16172836 : t42 = LOG(t41)
2303 16172836 : t43 = t42 + 0.2e1_dp
2304 16172836 : t44 = br_d1*t3
2305 16172836 : t45 = 0.1e1_dp/t19
2306 16172836 : t46 = t8*t45
2307 16172836 : t49 = br_d2*t24
2308 16172836 : t52 = pi**2
2309 16172836 : t53 = br_d3*t52
2310 16172836 : t54 = t26**2
2311 16172836 : t55 = t54*rho
2312 16172836 : t57 = 0.1e1_dp/t29/t19
2313 16172836 : t58 = t55*t57
2314 16172836 : t61 = t3*t52
2315 16172836 : t62 = br_d4*t61
2316 16172836 : t63 = t54*t26
2317 16172836 : t64 = t7*t63
2318 16172836 : t65 = t29**2
2319 16172836 : t66 = 0.1e1_dp/t65
2320 16172836 : t67 = t64*t66
2321 16172836 : t71 = t2*t52*pi
2322 16172836 : t72 = br_d5*t71
2323 16172836 : t73 = t54**2
2324 16172836 : t74 = t6*t73
2325 16172836 : t76 = 0.1e1_dp/t65/t19
2326 16172836 : t77 = t74*t76
2327 : t80 = br_d0 + 0.2e1_dp/0.3e1_dp*t44*t46 + 0.4e1_dp/0.9e1_dp*t49*t31 &
2328 : + 0.8e1_dp/0.27e2_dp*t53*t58 + 0.16e2_dp/0.81e2_dp*t62*t67 + 0.32e2_dp &
2329 16172836 : /0.243e3_dp*t72*t77
2330 16172836 : t81 = t43*t80
2331 16172836 : t82 = br_e1*t3
2332 16172836 : t85 = br_e2*t24
2333 16172836 : t88 = br_e3*t52
2334 16172836 : t91 = br_e4*t61
2335 16172836 : t94 = br_e5*t71
2336 : t97 = br_e0 + 0.2e1_dp/0.3e1_dp*t82*t46 + 0.4e1_dp/0.9e1_dp*t85*t31 &
2337 : + 0.8e1_dp/0.27e2_dp*t88*t58 + 0.16e2_dp/0.81e2_dp*t91*t67 + 0.32e2_dp &
2338 16172836 : /0.243e3_dp*t94*t77
2339 16172836 : t98 = 0.1e1_dp/t97
2340 16172836 : t99 = t81*t98
2341 16172836 : t100 = 8**(0.1e1_dp/0.3e1_dp)
2342 16172836 : t101 = t43**2
2343 16172836 : t102 = t101*t43
2344 16172836 : t103 = t80**2
2345 16172836 : t104 = t103*t80
2346 16172836 : t105 = t102*t104
2347 16172836 : t106 = t97**2
2348 16172836 : t108 = 0.1e1_dp/t106/t97
2349 16172836 : t109 = t105*t108
2350 16172836 : t110 = EXP(-t99)
2351 16172836 : t111 = 0.1e1_dp/0.3141592654e1_dp
2352 16172836 : t112 = t110*t111
2353 16172836 : t114 = t109*t112*t13
2354 16172836 : t115 = t114**(0.1e1_dp/0.3e1_dp)
2355 16172836 : t116 = 0.1e1_dp/t115
2356 16172836 : t117 = REAL(t100, KIND=dp)*t116
2357 16172836 : t118 = t117*R
2358 16172836 : t119 = t99*t118
2359 16172836 : t121 = EXP(0.2e1_dp*t119)
2360 16172836 : t123 = t121*R
2361 16172836 : t124 = t123*t43
2362 16172836 : t125 = t80*t98
2363 16172836 : t126 = t125*t117
2364 16172836 : t128 = t121*t43
2365 16172836 : t130 = t99 + t119
2366 16172836 : t131 = EXP(t130)
2367 : t133 = 0.2e1_dp*t121 - t124*t126 + t128*t125 + 0.2e1_dp + t99 + t119 &
2368 16172836 : - 0.4e1_dp*t131
2369 16172836 : t134 = rho*t133
2370 16172836 : t136 = EXP(-t130)
2371 16172836 : t138 = t136*REAL(t100, KIND=dp)*t116
2372 16172836 : e_0 = e_0 + (t134*t138/0.8e1_dp)*sx
2373 16172836 : IF (grad_deriv >= 1) THEN
2374 8506079 : t144 = 0.1e1_dp/t7/t26
2375 8506079 : t151 = 0.1e1_dp/t7/t27*gamma*t12
2376 8506079 : t154 = 0.1e1_dp/t35
2377 8506079 : t157 = t6*t26
2378 8506079 : t158 = t157*t30
2379 8506079 : t161 = t6*rho
2380 8506079 : t163 = t57*gamma
2381 8506079 : t164 = t163*t12
2382 8506079 : t175 = t36*t4
2383 : t178 = -0.2500000000e1_dp*t5*t144*t19 - 0.1250000000e0_dp*t5*t151 &
2384 : + 0.3e1_dp/0.4e1_dp*t154*t1*t4*t20*(0.40e2_dp/0.27e2_dp*t25 &
2385 : *t158 + 0.2e1_dp/0.27e2_dp*t25*t161*t164) - 0.5e1_dp/0.2e1_dp*t36 &
2386 8506079 : *t4*t144*t19 - t175*t151/0.8e1_dp
2387 8506079 : t179 = 0.1e1_dp/t41
2388 8506079 : t180 = t178*t179
2389 8506079 : t182 = t98*REAL(t100, KIND=dp)
2390 8506079 : t184 = t182*t116*R
2391 8506079 : t185 = t180*t80*t184
2392 8506079 : t187 = t7*t45
2393 8506079 : t190 = 0.1e1_dp/t6
2394 8506079 : t192 = t30*gamma
2395 8506079 : t193 = t192*t12
2396 8506079 : t201 = t54*t57
2397 8506079 : t205 = t66*gamma
2398 8506079 : t206 = t205*t12
2399 8506079 : t209 = t7*t55
2400 8506079 : t210 = t209*t66
2401 8506079 : t213 = t7*t54
2402 8506079 : t215 = t76*gamma
2403 8506079 : t216 = t215*t12
2404 8506079 : t220 = t6*t54*t27
2405 8506079 : t221 = t220*t76
2406 8506079 : t224 = t6*t63
2407 8506079 : t227 = 0.1e1_dp/t65/t29
2408 8506079 : t228 = t227*gamma
2409 8506079 : t229 = t228*t12
2410 : t232 = 0.10e2_dp/0.9e1_dp*t44*t187 + t44*t190*t193/0.18e2_dp + 0.40e2_dp &
2411 : /0.27e2_dp*t49*t158 + 0.2e1_dp/0.27e2_dp*t49*t161*t164 + &
2412 : 0.40e2_dp/0.27e2_dp*t53*t201 + 0.2e1_dp/0.27e2_dp*t53*t27*t206 &
2413 : + 0.320e3_dp/0.243e3_dp*t62*t210 + 0.16e2_dp/0.243e3_dp*t62*t213* &
2414 : t216 + 0.800e3_dp/0.729e3_dp*t72*t221 + 0.40e2_dp/0.729e3_dp*t72* &
2415 8506079 : t224*t229
2416 8506079 : t234 = t43*t232*t98
2417 8506079 : t235 = t234*t118
2418 8506079 : t237 = 0.1e1_dp/t106
2419 8506079 : t238 = t81*t237
2420 : t264 = 0.10e2_dp/0.9e1_dp*t82*t187 + t82*t190*t193/0.18e2_dp + 0.40e2_dp &
2421 : /0.27e2_dp*t85*t158 + 0.2e1_dp/0.27e2_dp*t85*t161*t164 + &
2422 : 0.40e2_dp/0.27e2_dp*t88*t201 + 0.2e1_dp/0.27e2_dp*t88*t27*t206 &
2423 : + 0.320e3_dp/0.243e3_dp*t91*t210 + 0.16e2_dp/0.243e3_dp*t91*t213* &
2424 : t216 + 0.800e3_dp/0.729e3_dp*t94*t221 + 0.40e2_dp/0.729e3_dp*t94* &
2425 8506079 : t224*t229
2426 8506079 : t267 = t238*t117*R*t264
2427 8506079 : t270 = 0.1e1_dp/t115/t114
2428 8506079 : t271 = REAL(t100, KIND=dp)*t270
2429 8506079 : t274 = t101*t104*t108*t110
2430 8506079 : t275 = t111*t13
2431 8506079 : t280 = t102*t103*t108
2432 8506079 : t285 = t106**2
2433 8506079 : t287 = t105/t285
2434 8506079 : t292 = t180*t125
2435 8506079 : t294 = t81*t237*t264
2436 : t302 = 0.3e1_dp*t274*t275*t180 + 0.3e1_dp*t280*t112*t13*t232 &
2437 : - 0.3e1_dp*t287*t112*t13*t264 + t109*(-t292 - t234 + t294) &
2438 8506079 : *t110*t275 - t109*t112/t26
2439 8506079 : t305 = t99*t271*R*t302
2440 : t308 = (0.2e1_dp*t185 + 0.2e1_dp*t235 - 0.2e1_dp*t267 - 0.2e1_dp/0.3e1_dp &
2441 8506079 : *t305)*t121
2442 8506079 : t310 = R*t43
2443 8506079 : t315 = t232*t98
2444 8506079 : t318 = t123*t81
2445 8506079 : t319 = t237*REAL(t100, KIND=dp)
2446 8506079 : t330 = t179*t80*t98
2447 8506079 : t333 = t80*t237
2448 8506079 : t336 = t305/0.3e1_dp
2449 8506079 : t337 = t292 + t234 - t294 + t185 + t235 - t267 - t336
2450 : t340 = 0.2e1_dp*t308 - t308*t310*t126 - t123*t180*t126 - t124 &
2451 : *t315*t117 + t318*t319*t116*t264 + t318*t182*t270* &
2452 : t302/0.3e1_dp + t308*t99 + t121*t178*t330 + t128*t315 - t128 &
2453 : *t333*t264 + t292 + t234 - t294 + t185 + t235 - t267 - t336 &
2454 8506079 : - 0.4e1_dp*t337*t131
2455 8506079 : t348 = t134*t136
2456 : e_rho = e_rho + (t133*t136*t117/0.8e1_dp + rho*t340*t138/0.8e1_dp - t134 &
2457 8506079 : *t337*t138/0.8e1_dp - t348*t271*t302/0.24e2_dp)*sx
2458 8506079 : t353 = t144*gamma*ndrho
2459 8506079 : t356 = t154*br_BB
2460 8506079 : t357 = t356*t3
2461 : t365 = 0.2500000000000000e0_dp*t5*t353 - t357*t7*t30*gamma*ndrho &
2462 8506079 : /0.9e1_dp + t175*t353/0.4e1_dp
2463 8506079 : t366 = t365*t179
2464 8506079 : t368 = t366*t80*t184
2465 8506079 : t371 = t192*ndrho
2466 8506079 : t375 = t163*ndrho
2467 8506079 : t379 = t205*ndrho
2468 8506079 : t383 = t215*ndrho
2469 8506079 : t387 = t228*ndrho
2470 : t390 = -t44*t7*t371/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t49*t157*t375 &
2471 : - 0.4e1_dp/0.27e2_dp*t53*t54*t379 - 0.32e2_dp/0.243e3_dp*t62*t209 &
2472 8506079 : *t383 - 0.80e2_dp/0.729e3_dp*t72*t220*t387
2473 8506079 : t392 = t43*t390*t98
2474 8506079 : t393 = t392*t118
2475 : t410 = -t82*t7*t371/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t85*t157*t375 &
2476 : - 0.4e1_dp/0.27e2_dp*t88*t54*t379 - 0.32e2_dp/0.243e3_dp*t91*t209 &
2477 8506079 : *t383 - 0.80e2_dp/0.729e3_dp*t94*t220*t387
2478 8506079 : t413 = t238*t117*R*t410
2479 8506079 : t426 = t366*t125
2480 8506079 : t428 = t81*t237*t410
2481 : t433 = 0.3e1_dp*t274*t275*t366 + 0.3e1_dp*t280*t112*t13*t390 &
2482 : - 0.3e1_dp*t287*t112*t13*t410 + t109*(-t426 - t392 + t428) &
2483 8506079 : *t110*t275
2484 8506079 : t436 = t99*t271*R*t433
2485 : t439 = (0.2e1_dp*t368 + 0.2e1_dp*t393 - 0.2e1_dp*t413 - 0.2e1_dp/0.3e1_dp &
2486 8506079 : *t436)*t121
2487 8506079 : t445 = t390*t98
2488 8506079 : t461 = t436/0.3e1_dp
2489 8506079 : t462 = t426 + t392 - t428 + t368 + t393 - t413 - t461
2490 : t465 = 0.2e1_dp*t439 - t439*t310*t126 - t123*t366*t126 - t124 &
2491 : *t445*t117 + t318*t319*t116*t410 + t318*t182*t270* &
2492 : t433/0.3e1_dp + t439*t99 + t121*t365*t330 + t128*t445 - t128 &
2493 : *t333*t410 + t426 + t392 - t428 + t368 + t393 - t413 - t461 &
2494 8506079 : - 0.4e1_dp*t462*t131
2495 : e_ndrho = e_ndrho + (rho*t465*t138/0.8e1_dp - t134*t462*t138/0.8e1_dp - t348 &
2496 8506079 : *t271*t433/0.24e2_dp)*sx
2497 8506079 : t479 = t8*t30
2498 8506079 : t480 = t479*gamma
2499 : t485 = -0.1000000000e1_dp*t5*t9*gamma + 0.4e1_dp/0.9e1_dp*t357*t480 &
2500 8506079 : - t36*t37*gamma
2501 8506079 : t486 = t485*t179
2502 8506079 : t488 = t486*t80*t184
2503 8506079 : t492 = t28*t57
2504 8506079 : t493 = t492*gamma
2505 8506079 : t496 = t55*t66
2506 8506079 : t497 = t496*gamma
2507 8506079 : t500 = t64*t76
2508 8506079 : t501 = t500*gamma
2509 8506079 : t504 = t74*t227
2510 8506079 : t505 = t504*gamma
2511 : t508 = 0.4e1_dp/0.9e1_dp*t44*t480 + 0.16e2_dp/0.27e2_dp*t49*t493 + &
2512 : 0.16e2_dp/0.27e2_dp*t53*t497 + 0.128e3_dp/0.243e3_dp*t62*t501 + 0.320e3_dp &
2513 8506079 : /0.729e3_dp*t72*t505
2514 8506079 : t510 = t43*t508*t98
2515 8506079 : t511 = t510*t118
2516 : t523 = 0.4e1_dp/0.9e1_dp*t82*t480 + 0.16e2_dp/0.27e2_dp*t85*t493 + &
2517 : 0.16e2_dp/0.27e2_dp*t88*t497 + 0.128e3_dp/0.243e3_dp*t91*t501 + 0.320e3_dp &
2518 8506079 : /0.729e3_dp*t94*t505
2519 8506079 : t526 = t238*t117*R*t523
2520 8506079 : t539 = t486*t125
2521 8506079 : t541 = t81*t237*t523
2522 : t546 = 0.3e1_dp*t274*t275*t486 + 0.3e1_dp*t280*t112*t13*t508 &
2523 : - 0.3e1_dp*t287*t112*t13*t523 + t109*(-t539 - t510 + t541) &
2524 8506079 : *t110*t275
2525 8506079 : t549 = t99*t271*R*t546
2526 : t552 = (0.2e1_dp*t488 + 0.2e1_dp*t511 - 0.2e1_dp*t526 - 0.2e1_dp/0.3e1_dp &
2527 8506079 : *t549)*t121
2528 8506079 : t558 = t508*t98
2529 8506079 : t574 = t549/0.3e1_dp
2530 8506079 : t575 = t539 + t510 - t541 + t488 + t511 - t526 - t574
2531 : t578 = 0.2e1_dp*t552 - t552*t310*t126 - t123*t486*t126 - t124 &
2532 : *t558*t117 + t318*t319*t116*t523 + t318*t182*t270* &
2533 : t546/0.3e1_dp + t552*t99 + t121*t485*t330 + t128*t558 - t128 &
2534 : *t333*t523 + t539 + t510 - t541 + t488 + t511 - t526 - t574 &
2535 8506079 : - 0.4e1_dp*t575*t131
2536 : e_tau = e_tau + (rho*t578*t138/0.8e1_dp - t134*t575*t138/0.8e1_dp - t348 &
2537 8506079 : *t271*t546/0.24e2_dp)*sx
2538 : t597 = 0.2500000000000000e0_dp*t5*t9 - t356*t3*t8*t30/0.9e1_dp &
2539 8506079 : + t36*t37/0.4e1_dp
2540 8506079 : t598 = t597*t179
2541 8506079 : t600 = t598*t80*t184
2542 : t612 = -t44*t479/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t49*t492 - 0.4e1_dp/ &
2543 : 0.27e2_dp*t53*t496 - 0.32e2_dp/0.243e3_dp*t62*t500 - 0.80e2_dp/0.729e3_dp &
2544 8506079 : *t72*t504
2545 8506079 : t614 = t43*t612*t98
2546 8506079 : t615 = t614*t118
2547 : t627 = -t82*t479/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t85*t492 - 0.4e1_dp/ &
2548 : 0.27e2_dp*t88*t496 - 0.32e2_dp/0.243e3_dp*t91*t500 - 0.80e2_dp/0.729e3_dp &
2549 8506079 : *t94*t504
2550 8506079 : t630 = t238*t117*R*t627
2551 8506079 : t643 = t598*t125
2552 8506079 : t645 = t81*t237*t627
2553 : t650 = 0.3e1_dp*t274*t275*t598 + 0.3e1_dp*t280*t112*t13*t612 &
2554 : - 0.3e1_dp*t287*t112*t13*t627 + t109*(-t643 - t614 + t645) &
2555 8506079 : *t110*t275
2556 8506079 : t653 = t99*t271*R*t650
2557 : t656 = (0.2e1_dp*t600 + 0.2e1_dp*t615 - 0.2e1_dp*t630 - 0.2e1_dp/0.3e1_dp &
2558 8506079 : *t653)*t121
2559 8506079 : t662 = t612*t98
2560 8506079 : t678 = t653/0.3e1_dp
2561 8506079 : t679 = t643 + t614 - t645 + t600 + t615 - t630 - t678
2562 : t682 = 0.2e1_dp*t656 - t656*t310*t126 - t123*t598*t126 - t124 &
2563 : *t662*t117 + t318*t319*t116*t627 + t318*t182*t270* &
2564 : t650/0.3e1_dp + t656*t99 + t121*t597*t330 + t128*t662 - t128 &
2565 : *t333*t627 + t643 + t614 - t645 + t600 + t615 - t630 - t678 &
2566 8506079 : - 0.4e1_dp*t679*t131
2567 : e_laplace_rho = e_laplace_rho + (rho*t682*t138/0.8e1_dp - t134*t679*t138/0.8e1_dp - t348 &
2568 8506079 : *t271*t650/0.24e2_dp)*sx
2569 : END IF
2570 :
2571 16172836 : END SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_lte_b
2572 :
2573 : END MODULE xc_xbecke_roussel
|
