(git:08825b9)
Loading...
Searching...
No Matches
xc_xbecke_roussel.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculates the 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! **************************************************************************************************
18MODULE 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
77CONTAINS
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 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 CALL cite_reference(beckeroussel1989)
96 CALL cite_reference(proynov2007)
97 IF (PRESENT(reference)) THEN
98 reference = "A.D. Becke, M.R. Roussel, "// &
99 "Phys. Rev. A, vol. 39, n 8, pp. 3761-3767, (1989) {spin unpolarized}"
100 END IF
101 IF (PRESENT(shortform)) THEN
102 shortform = "Becke-Roussel exchange hole (spin unpolarized)"
103 END IF
104 IF (PRESENT(needs)) THEN
105 needs%rho = .true.
106 needs%norm_drho = .true.
107 needs%tau = .true.
108 needs%laplace_rho = .true.
109 END IF
110
111 IF (PRESENT(max_deriv)) max_deriv = 1
112
113 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 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 CALL cite_reference(beckeroussel1989)
132 CALL cite_reference(proynov2007)
133 IF (PRESENT(reference)) THEN
134 reference = "A.D. Becke, M.R. Roussel, "// &
135 "Phys. Rev. A, vol. 39, n 8, pp. 3761-3767, (1989) {spin polarized}"
136 END IF
137 IF (PRESENT(shortform)) THEN
138 shortform = "Becke-Roussel exchange hole (spin polarized)"
139 END IF
140 IF (PRESENT(needs)) THEN
141 needs%rho_spin = .true.
142 needs%norm_drho_spin = .true.
143 needs%tau_spin = .true.
144 needs%laplace_rho_spin = .true.
145 END IF
146 IF (PRESENT(max_deriv)) max_deriv = 1
147
148 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 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 POINTER :: dummy, e_0, e_laplace_rho, e_ndrho, &
175 e_rho, e_tau, laplace_rho, norm_drho, &
176 rho, tau
177 TYPE(xc_derivative_type), POINTER :: deriv
178
179 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 rho_cutoff=epsilon_rho)
184 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
185
186 dummy => rho
187
188 e_0 => dummy
189 e_rho => dummy
190 e_ndrho => dummy
191 e_tau => dummy
192 e_laplace_rho => dummy
193
194 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
195 allocate_deriv=.true.)
196 CALL xc_derivative_get(deriv, deriv_data=e_0)
197 IF (grad_deriv >= 1) THEN
198 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
199 allocate_deriv=.true.)
200 CALL xc_derivative_get(deriv, deriv_data=e_rho)
201 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
202 allocate_deriv=.true.)
203 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
204 deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
205 allocate_deriv=.true.)
206 CALL xc_derivative_get(deriv, deriv_data=e_tau)
207 deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho], &
208 allocate_deriv=.true.)
209 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho)
210 END IF
211 IF (grad_deriv > 1) THEN
212 cpabort("derivatives bigger than 1 not implemented")
213 END IF
214
215 CALL section_vals_val_get(br_params, "scale_x", r_val=sx)
216 CALL section_vals_val_get(br_params, "CUTOFF_RADIUS", r_val=r)
217 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!$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 CALL timestop(handle)
234 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 SUBROUTINE xbecke_roussel_lda_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
262 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 my_rho = 0.5_dp*max(rho(ip), 0.0_dp)
281 IF (my_rho > epsilon_rho) THEN
282 my_ndrho = 0.5_dp*max(norm_drho(ip), epsilon(0.0_dp)*1.e4_dp)
283 my_tau = 0.5_dp*max(epsilon(0.0_dp)*1.e4_dp, tau(ip))
284 my_laplace_rho = 0.5_dp*laplace_rho(ip)
285
286 t1 = pi**(0.1e1_dp/0.3e1_dp)
287 t2 = t1**2
288 t3 = my_rho**(0.1e1_dp/0.3e1_dp)
289 t4 = t3**2
290 t5 = t4*my_rho
291 t8 = my_ndrho**2
292 t9 = 0.1e1_dp/my_rho
293 ! *** CP2K defines tau in a different way as compared to Becke !!!
294 t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
295 t16 = 0.1e1_dp/t15
296 yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
297 IF (r == 0.0_dp) THEN
298 IF (yval <= 0.0_dp) THEN
299 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 sx, gamma, grad_deriv)
303 ! VERY UGLY HACK e_0 has to multiplied by the factor 2
304 e_diff = e_0(ip) - e_old
305 e_0(ip) = e_0(ip) + e_diff
306 ELSE
307 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 sx, gamma, grad_deriv)
311 ! VERY UGLY HACK e_0 has to multiplied by the factor 2
312 e_diff = e_0(ip) - e_old
313 e_0(ip) = e_0(ip) + e_diff
314 END IF
315 ELSE
316 IF (yval <= 0.0_dp) THEN
317 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 sx, r, gamma, grad_deriv)
321 ! VERY UGLY HACK e_0 has to multiplied by the factor 2
322 e_diff = e_0(ip) - e_old
323 e_0(ip) = e_0(ip) + e_diff
324 ELSE
325 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 sx, r, gamma, grad_deriv)
329 ! VERY UGLY HACK e_0 has to multiplied by the factor 2
330 e_diff = e_0(ip) - e_old
331 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 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 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 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, &
365 e_laplace_rhob, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, e_tau_a, e_tau_b, laplace_rhoa, &
366 laplace_rhob, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
367 TYPE(xc_derivative_type), POINTER :: deriv
368
369 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 rho_cutoff=epsilon_rho)
375 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
376
377 dummy => rhoa
378
379 e_0 => dummy
380 e_rhoa => dummy
381 e_rhob => dummy
382 e_ndrhoa => dummy
383 e_ndrhob => dummy
384 e_tau_a => dummy
385 e_tau_b => dummy
386 e_laplace_rhoa => dummy
387 e_laplace_rhob => dummy
388
389 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
390 allocate_deriv=.true.)
391 CALL xc_derivative_get(deriv, deriv_data=e_0)
392 IF (grad_deriv >= 1) THEN
393 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
394 allocate_deriv=.true.)
395 CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
396 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
397 allocate_deriv=.true.)
398 CALL xc_derivative_get(deriv, deriv_data=e_rhob)
399
400 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
401 allocate_deriv=.true.)
402 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
403 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
404 allocate_deriv=.true.)
405 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
406
407 deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], &
408 allocate_deriv=.true.)
409 CALL xc_derivative_get(deriv, deriv_data=e_tau_a)
410 deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], &
411 allocate_deriv=.true.)
412 CALL xc_derivative_get(deriv, deriv_data=e_tau_b)
413
414 deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa], &
415 allocate_deriv=.true.)
416 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa)
417 deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob], &
418 allocate_deriv=.true.)
419 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob)
420 END IF
421 IF (grad_deriv > 1) THEN
422 cpabort("derivatives bigger than 1 not implemented")
423 END IF
424
425 CALL section_vals_val_get(br_params, "scale_x", r_val=sx)
426 CALL section_vals_val_get(br_params, "CUTOFF_RADIUS", r_val=r)
427 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!$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 CALL timestop(handle)
452 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 SUBROUTINE xbecke_roussel_lsd_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
478 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 my_rho = max(rho(ip), 0.0_dp)
497 IF (my_rho > epsilon_rho) THEN
498 my_ndrho = max(norm_drho(ip), epsilon(0.0_dp)*1.e4_dp)
499 my_tau = 1.0_dp*max(epsilon(0.0_dp)*1.e4_dp, tau(ip))
500 my_laplace_rho = 1.0_dp*laplace_rho(ip)
501
502 t1 = pi**(0.1e1_dp/0.3e1_dp)
503 t2 = t1**2
504 t3 = my_rho**(0.1e1_dp/0.3e1_dp)
505 t4 = t3**2
506 t5 = t4*my_rho
507 t8 = my_ndrho**2
508 t9 = 0.1e1_dp/my_rho
509 ! *** CP2K defines tau in a different way as compared to Becke !!!
510 t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
511 t16 = 0.1e1_dp/t15
512 yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
513 IF (r == 0.0_dp) THEN
514 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 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 sx, gamma, grad_deriv)
522 END IF
523 ELSE
524 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 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 sx, r, gamma, grad_deriv)
532 END IF
533 END IF
534 END IF
535 END DO
536
537!$OMP END DO
538
539 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 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 t1 = pi**(0.1e1_dp/0.3e1_dp)
580 t2 = t1**2
581 t3 = rho**(0.1e1_dp/0.3e1_dp)
582 t4 = t3**2
583 t5 = t4*rho
584 t9 = ndrho**2
585 t16 = laplace_rho/0.6e1_dp - gamma*(real(2*tau, kind=dp) - t9/rho/0.4e1_dp)/0.3e1_dp
586 t17 = 0.1e1_dp/t16
587 yval = 0.2e1_dp/0.3e1_dp*t2*t5*t17
588 t19 = t3*rho
589 t20 = 0.3141592654e1_dp**(0.1e1_dp/0.3e1_dp)
590 t21 = t19*t20
591 t22 = br_a1*t2
592 t23 = t5*t17
593 t26 = 0.2e1_dp/0.3e1_dp*t22*t23 + br_a2
594 t27 = atan(t26)
595 t28 = -t27 + br_a3
596 t29 = br_c1*t2
597 t32 = t1*pi
598 t33 = br_c2*t32
599 t34 = rho**2
600 t35 = t34*rho
601 t36 = t3*t35
602 t37 = t16**2
603 t38 = 0.1e1_dp/t37
604 t39 = t36*t38
605 t42 = pi**2
606 t43 = br_c3*t42
607 t44 = t34**2
608 t45 = t44*rho
609 t47 = 0.1e1_dp/t37/t16
610 t48 = t45*t47
611 t51 = t2*t42
612 t52 = br_c4*t51
613 t53 = t44*t34
614 t54 = t4*t53
615 t55 = t37**2
616 t56 = 0.1e1_dp/t55
617 t57 = t54*t56
618 t61 = t1*t42*pi
619 t62 = br_c5*t61
620 t63 = t44**2
621 t64 = t3*t63
622 t66 = 0.1e1_dp/t55/t16
623 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 /0.243e3_dp*t62*t67
627 t71 = t28*t70
628 t72 = br_b1*t2
629 t75 = br_b2*t32
630 t78 = br_b3*t42
631 t81 = br_b4*t51
632 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 /0.243e3_dp*t84*t67
636 t88 = 0.1e1_dp/t87
637 t89 = t71*t88
638 t91 = exp(t89/0.3e1_dp)
639 t92 = t21*t91
640 t93 = 0.1e1_dp/t28
641 t94 = 0.1e1_dp/t70
642 t95 = t93*t94
643 t96 = exp(-t89)
644 t97 = t88*t96
645 t100 = 0.1e1_dp - t96 - t71*t97/0.2e1_dp
646 t101 = t87*t100
647 t102 = t95*t101
648 e_0 = e_0 + (-t92*t102)*sx
649 IF (grad_deriv >= 1) THEN
650 t108 = t4*t17
651 t111 = 0.1e1_dp/t3
652 t113 = t38*gamma
653 t114 = t113*t9
654 t117 = 0.10e2_dp/0.9e1_dp*t22*t108 + t22*t111*t114/0.18e2_dp
655 t118 = t26**2
656 t120 = 0.1e1_dp/(0.1e1_dp + t118)
657 t121 = t117*t120
658 t122 = t70*t88
659 t129 = t3*t34
660 t130 = t129*t38
661 t134 = t47*gamma
662 t135 = t134*t9
663 t138 = t44*t47
664 t142 = t56*gamma
665 t143 = t142*t9
666 t146 = t4*t45
667 t147 = t146*t56
668 t150 = t4*t44
669 t152 = t66*gamma
670 t153 = t152*t9
671 t157 = t3*t44*t35
672 t158 = t157*t66
673 t161 = t3*t53
674 t164 = 0.1e1_dp/t55/t37
675 t165 = t164*gamma
676 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 *t166
683 t170 = t28*t169
684 t172 = t87**2
685 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 *t166
692 t202 = -t121*t122 + t170*t88 - t71*t173*t199
693 t207 = t28**2
694 t208 = 0.1e1_dp/t207
695 t209 = t91*t208
696 t217 = t21*t91*t93
697 t218 = t70**2
698 t220 = 0.1e1_dp/t218*t87
699 t227 = -t202
700 t229 = t122*t96
701 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 *t199/0.2e1_dp - t71*t88*t227*t96/0.2e1_dp))*sx
707 t246 = t4*t38
708 t249 = t120*t70
709 t252 = t22*t246*gamma*ndrho*t249*t88
710 t255 = t113*ndrho
711 t259 = t134*ndrho
712 t263 = t142*ndrho
713 t267 = t152*ndrho
714 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 *t267 - 0.80e2_dp/0.729e3_dp*t62*t157*t271
718 t275 = t28*t274
719 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 *t267 - 0.80e2_dp/0.729e3_dp*t84*t157*t271
723 t295 = t71*t173*t293
724 t304 = t208*t94*t87
725 t307 = t100*br_a1*t2
726 t308 = ndrho*t120
727 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 /0.2e1_dp))*sx
734 t340 = t5*t38
735 t341 = t22*t340
736 t342 = gamma*t120
737 t344 = t341*t342*t122
738 t346 = t340*gamma
739 t349 = t36*t47
740 t350 = t349*gamma
741 t353 = t45*t56
742 t354 = t353*gamma
743 t357 = t54*t66
744 t358 = t357*gamma
745 t361 = t64*t164
746 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 /0.729e3_dp*t62*t362
750 t366 = t28*t365
751 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 /0.729e3_dp*t84*t362
755 t381 = t71*t173*t379
756 t387 = t35*t20
757 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 *t96/0.2e1_dp))*sx
764 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 *t62*t361
768 t435 = t28*t434
769 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 *t84*t361
773 t450 = t71*t173*t448
774 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 + t71*t234*t448/0.2e1_dp - t71*t88*t471*t96/0.2e1_dp))*sx
780 END IF
781 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 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 t1 = pi**(0.1e1_dp/0.3e1_dp)
822 t2 = t1**2
823 t3 = rho**(0.1e1_dp/0.3e1_dp)
824 t4 = t3**2
825 t5 = t4*rho
826 t6 = t2*t5
827 t9 = ndrho**2
828 t16 = laplace_rho/0.6e1_dp - gamma*(real(2*tau, kind=dp) - t9/rho/0.4e1_dp)/0.3e1_dp
829 t17 = 0.1e1_dp/t16
830 yval = 0.2e1_dp/0.3e1_dp*t6*t17
831 t19 = t3*rho
832 t20 = 0.3141592654e1_dp**(0.1e1_dp/0.3e1_dp)
833 t21 = t19*t20
834 t22 = 0.1e1_dp/br_bb
835 t23 = 0.1e1_dp/t2
836 t24 = t22*t23
837 t25 = 0.1e1_dp/t5
838 t26 = t25*t16
839 t29 = br_bb**2
840 t30 = t1*pi
841 t31 = t29*t30
842 t32 = rho**2
843 t33 = t32*rho
844 t34 = t3*t33
845 t35 = t16**2
846 t36 = 0.1e1_dp/t35
847 t37 = t34*t36
848 t41 = sqrt(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t31*t37)
849 t42 = t41*t22
850 t43 = t23*t25
851 t47 = 0.1500000000000000e1_dp*t24*t26 + 0.3e1_dp/0.2e1_dp*t42*t43 &
852 *t16
853 t48 = log(t47)
854 t49 = t48 + 0.2e1_dp
855 t50 = br_d1*t2
856 t51 = t5*t17
857 t54 = br_d2*t30
858 t57 = pi**2
859 t58 = br_d3*t57
860 t59 = t32**2
861 t60 = t59*rho
862 t62 = 0.1e1_dp/t35/t16
863 t63 = t60*t62
864 t66 = t2*t57
865 t67 = br_d4*t66
866 t68 = t59*t32
867 t69 = t4*t68
868 t70 = t35**2
869 t71 = 0.1e1_dp/t70
870 t72 = t69*t71
871 t76 = t1*t57*pi
872 t77 = br_d5*t76
873 t78 = t59**2
874 t79 = t3*t78
875 t81 = 0.1e1_dp/t70/t16
876 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 /0.243e3_dp*t77*t82
880 t86 = t49*t85
881 t87 = br_e1*t2
882 t90 = br_e2*t30
883 t93 = br_e3*t57
884 t96 = br_e4*t66
885 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 /0.243e3_dp*t99*t82
889 t103 = 0.1e1_dp/t102
890 t104 = t86*t103
891 t106 = exp(t104/0.3e1_dp)
892 t107 = t21*t106
893 t108 = 0.1e1_dp/t49
894 t109 = 0.1e1_dp/t85
895 t110 = t108*t109
896 t111 = exp(-t104)
897 t112 = t103*t111
898 t115 = 0.1e1_dp - t111 - t86*t112/0.2e1_dp
899 t117 = t110*t102*t115
900 e_0 = e_0 + (-t107*t117)*sx
901 IF (grad_deriv >= 1) THEN
902 t124 = 0.1e1_dp/t4/t32
903 t131 = 0.1e1_dp/t4/t33*gamma*t9
904 t134 = 0.1e1_dp/t41
905 t137 = t3*t32
906 t138 = t137*t36
907 t142 = t62*gamma
908 t143 = t142*t9
909 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 t42*t23*t124*t16 - t154*t131/0.8e1_dp
914 t158 = 0.1e1_dp/t47
915 t159 = t157*t158
916 t160 = t85*t103
917 t162 = t4*t17
918 t165 = 0.1e1_dp/t3
919 t167 = t36*gamma
920 t168 = t167*t9
921 t176 = t59*t62
922 t180 = t71*gamma
923 t181 = t180*t9
924 t184 = t4*t60
925 t185 = t184*t71
926 t188 = t4*t59
927 t190 = t81*gamma
928 t191 = t190*t9
929 t195 = t3*t59*t33
930 t196 = t195*t81
931 t199 = t3*t68
932 t202 = 0.1e1_dp/t70/t35
933 t203 = t202*gamma
934 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 *t204
941 t208 = t49*t207
942 t210 = t102**2
943 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 *t204
950 t240 = t159*t160 + t208*t103 - t86*t211*t237
951 t245 = t49**2
952 t248 = t21*t106/t245
953 t249 = t109*t102
954 t255 = t21*t106*t108
955 t256 = t85**2
956 t258 = 0.1e1_dp/t256*t102
957 t265 = -t240
958 t267 = t160*t111
959 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 *t237/0.2e1_dp - t86*t103*t265*t111/0.2e1_dp))*sx
965 t285 = t124*gamma*ndrho
966 t288 = t134*br_bb
967 t289 = t288*t2
968 t297 = 0.2500000000000000e0_dp*t24*t285 - t289*t4*t36*gamma* &
969 ndrho/0.9e1_dp + t154*t285/0.4e1_dp
970 t298 = t297*t158
971 t301 = t167*ndrho
972 t305 = t142*ndrho
973 t309 = t180*ndrho
974 t313 = t190*ndrho
975 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 *t313 - 0.80e2_dp/0.729e3_dp*t77*t195*t317
979 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 *t313 - 0.80e2_dp/0.729e3_dp*t99*t195*t317
983 t341 = t298*t160 + t321*t103 - t86*t211*t338
984 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 t111/0.2e1_dp))*sx
990 t376 = t5*t36
991 t377 = t376*gamma
992 t382 = -0.1000000000e1_dp*t24*t25*gamma + 0.4e1_dp/0.9e1_dp*t289*t377 &
993 - t42*t43*gamma
994 t383 = t382*t158
995 t387 = t34*t62
996 t388 = t387*gamma
997 t391 = t60*t71
998 t392 = t391*gamma
999 t395 = t69*t81
1000 t396 = t395*gamma
1001 t399 = t79*t202
1002 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 /0.729e3_dp*t77*t400
1006 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 /0.729e3_dp*t99*t400
1010 t419 = t383*t160 + t404*t103 - t86*t211*t416
1011 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 /0.2e1_dp))*sx
1017 t458 = 0.2500000000000000e0_dp*t24*t25 - t288*t6*t36/0.9e1_dp + &
1018 t42*t43/0.4e1_dp
1019 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 *t77*t399
1023 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 *t99*t399
1027 t487 = t459*t160 + t472*t103 - t86*t211*t484
1028 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 /0.2e1_dp))*sx
1034 END IF
1035
1036 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 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 t1 = 8**(0.1e1_dp/0.3e1_dp)
1072 t2 = t1**2
1073 t3 = pi**(0.1e1_dp/0.3e1_dp)
1074 t4 = t3**2
1075 t6 = rho**(0.1e1_dp/0.3e1_dp)
1076 t7 = t6**2
1077 t11 = ndrho**2
1078 t12 = 0.1e1_dp/rho
1079 t18 = laplace_rho/0.6e1_dp - gamma*(real(2*tau, kind=dp) - t11*t12/0.4e1_dp)/0.3e1_dp
1080 t20 = t7*rho/t18
1081 t24 = atan(0.2e1_dp/0.3e1_dp*br_a1*t4*t20 + br_a2)
1082 t25 = -t24 + br_a3
1083 t26 = t25**2
1084 t31 = t3*pi
1085 t33 = rho**2
1086 t36 = t18**2
1087 t38 = t6*t33*rho/t36
1088 t41 = pi**2
1089 t43 = t33**2
1090 t47 = t43*rho/t36/t18
1091 t50 = t4*t41
1092 t54 = t36**2
1093 t56 = t7*t43*t33/t54
1094 t60 = t3*t41*pi
1095 t62 = t43**2
1096 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 *br_c4*t50*t56 + 0.32e2_dp/0.243e3_dp*br_c5*t60*t66
1100 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 *br_b4*t50*t56 + 0.32e2_dp/0.243e3_dp*br_b5*t60*t66
1104 t89 = t88**2
1105 t96 = exp(-t25*t69/t88)
1106 t101 = (t26*t25*t70*t69/t89/t88*t96/0.3141592654e1_dp* &
1107 t12)**(0.1e1_dp/0.3e1_dp)
1108 brval = real(t2, kind=dp)*t101/0.8e1_dp
1109
1110 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 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 sx, r, gamma, grad_deriv)
1118 END IF
1119
1120 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 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 t1 = pi**(0.1e1_dp/0.3e1_dp)
1167 t2 = t1**2
1168 t3 = br_a1*t2
1169 t4 = rho**(0.1e1_dp/0.3e1_dp)
1170 t5 = t4**2
1171 t6 = t5*rho
1172 t9 = ndrho**2
1173 t10 = 0.1e1_dp/rho
1174 t16 = laplace_rho/0.6e1_dp - gamma*(real(2*tau, kind=dp) - t9*t10/0.4e1_dp)/0.3e1_dp
1175 t17 = 0.1e1_dp/t16
1176 t18 = t6*t17
1177 t21 = 0.2e1_dp/0.3e1_dp*t3*t18 + br_a2
1178 t22 = atan(t21)
1179 t23 = -t22 + br_a3
1180 t24 = br_c1*t2
1181 t27 = t1*pi
1182 t28 = br_c2*t27
1183 t29 = rho**2
1184 t30 = t29*rho
1185 t31 = t4*t30
1186 t32 = t16**2
1187 t33 = 0.1e1_dp/t32
1188 t34 = t31*t33
1189 t37 = pi**2
1190 t38 = br_c3*t37
1191 t39 = t29**2
1192 t40 = t39*rho
1193 t42 = 0.1e1_dp/t32/t16
1194 t43 = t40*t42
1195 t46 = t2*t37
1196 t47 = br_c4*t46
1197 t48 = t39*t29
1198 t49 = t5*t48
1199 t50 = t32**2
1200 t51 = 0.1e1_dp/t50
1201 t52 = t49*t51
1202 t56 = t1*t37*pi
1203 t57 = br_c5*t56
1204 t58 = t39**2
1205 t59 = t4*t58
1206 t61 = 0.1e1_dp/t50/t16
1207 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 /0.243e3_dp*t57*t62
1211 t66 = t23*t65
1212 t67 = br_b1*t2
1213 t70 = br_b2*t27
1214 t73 = br_b3*t37
1215 t76 = br_b4*t46
1216 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 /0.243e3_dp*t79*t62
1220 t83 = 0.1e1_dp/t82
1221 t84 = t66*t83
1222 t85 = 8**(0.1e1_dp/0.3e1_dp)
1223 t86 = t23**2
1224 t87 = t86*t23
1225 t88 = t65**2
1226 t89 = t88*t65
1227 t90 = t87*t89
1228 t91 = t82**2
1229 t93 = 0.1e1_dp/t91/t82
1230 t94 = t90*t93
1231 t95 = exp(-t84)
1232 t96 = 0.1e1_dp/0.3141592654e1_dp
1233 t97 = t95*t96
1234 t99 = t94*t97*t10
1235 t100 = t99**(0.1e1_dp/0.3e1_dp)
1236 t101 = 0.1e1_dp/t100
1237 t102 = real(t85, kind=dp)*t101
1238 t103 = t102*r
1239 t104 = t84*t103
1240 t106 = exp(t84 - t104)
1241 t108 = t106*t23
1242 t109 = t65*t83
1243 t110 = t108*t109
1244 t114 = t83*real(t85, kind=dp)*t101*r
1245 t117 = exp(-t84 - t104)
1246 t119 = t117*t23
1247 t120 = t119*t109
1248 t123 = -0.2e1_dp*t106 + t110 - t108*t65*t114 + 0.2e1_dp*t117 + t120 &
1249 + t119*t65*t114
1250 t124 = rho*t123
1251 e_0 = e_0 + (t124*t102/0.8e1_dp)*sx
1252 IF (grad_deriv >= 1) THEN
1253 t129 = t5*t17
1254 t132 = 0.1e1_dp/t4
1255 t134 = t33*gamma
1256 t135 = t134*t9
1257 t138 = 0.10e2_dp/0.9e1_dp*t3*t129 + t3*t132*t135/0.18e2_dp
1258 t139 = t21**2
1259 t141 = 0.1e1_dp/(0.1e1_dp + t139)
1260 t142 = t138*t141
1261 t143 = t142*t109
1262 t149 = t4*t29
1263 t150 = t149*t33
1264 t153 = t4*rho
1265 t155 = t42*gamma
1266 t156 = t155*t9
1267 t159 = t39*t42
1268 t163 = t51*gamma
1269 t164 = t163*t9
1270 t167 = t5*t40
1271 t168 = t167*t51
1272 t171 = t5*t39
1273 t173 = t61*gamma
1274 t174 = t173*t9
1275 t178 = t4*t39*t30
1276 t179 = t178*t61
1277 t182 = t4*t48
1278 t185 = 0.1e1_dp/t50/t32
1279 t186 = t185*gamma
1280 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 t182*t187
1287 t192 = t23*t190*t83
1288 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 t182*t187
1295 t221 = t66*t193*t219
1296 t223 = t142*t65*t114
1297 t224 = t192*t103
1298 t225 = t66*t193
1299 t227 = t102*r*t219
1300 t228 = t225*t227
1301 t231 = real(t85, kind=dp)/t100/t99
1302 t232 = t86*t89
1303 t233 = t93*t95
1304 t235 = t96*t10
1305 t240 = t87*t88*t93
1306 t245 = t91**2
1307 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 t221)*t95*t235 - t94*t97/t29
1311 t261 = t231*r*t259
1312 t263 = t84*t261/0.3e1_dp
1313 t265 = (-t143 + t192 - t221 + t223 - t224 + t228 + t263)*t106
1314 t268 = t106*t138
1315 t269 = t141*t65
1316 t270 = t269*t83
1317 t272 = t190*t83
1318 t274 = t65*t193
1319 t275 = t274*t219
1320 t283 = t108*t274
1321 t288 = (t143 - t192 + t221 + t223 - t224 + t228 + t263)*t117
1322 t291 = t117*t138
1323 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 /0.3e1_dp
1330 e_rho = e_rho + (t123*real(t85, kind=dp)*t101/0.8e1_dp + rho*t305*t102/0.8e1_dp &
1331 - t124*t231*t259/0.24e2_dp)*sx
1332 t312 = t5*t33
1333 t314 = gamma*ndrho
1334 t315 = t314*t270
1335 t317 = t3*t312*t315/0.9e1_dp
1336 t319 = t134*ndrho
1337 t323 = t155*ndrho
1338 t327 = t163*ndrho
1339 t331 = t173*ndrho
1340 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 *t331 - 0.80e2_dp/0.729e3_dp*t57*t178*t335
1344 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 *t331 - 0.80e2_dp/0.729e3_dp*t79*t178*t335
1348 t358 = t66*t193*t356
1349 t359 = t3*t5
1350 t361 = t270*t103
1351 t363 = t359*t319*t361/0.9e1_dp
1352 t364 = t340*t103
1353 t366 = t102*r*t356
1354 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 t94*(-t317 - t340 + t358)*t95*t235
1358 t390 = t231*r*t388
1359 t392 = t84*t390/0.3e1_dp
1360 t394 = (t317 + t340 - t358 - t363 - t364 + t367 + t392)*t106
1361 t397 = t106*br_a1
1362 t399 = t2*t5*t33
1363 t403 = t338*t83
1364 t405 = t274*t356
1365 t409 = t397*t2
1366 t410 = t312*gamma
1367 t414 = ndrho*t141*t65*t114
1368 t423 = (-t317 - t340 + t358 - t363 - t364 + t367 + t392)*t117
1369 t426 = t117*br_a1
1370 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 + t119*t338*t114 - t301*t366 - t120*t390/0.3e1_dp
1377 e_ndrho = e_ndrho + (rho*t443*t102/0.8e1_dp - t124*t231*t388/0.24e2_dp)*sx
1378 t450 = t6*t33
1379 t455 = 0.4e1_dp/0.9e1_dp*t3*t450*gamma*t141*t109
1380 t456 = t450*gamma
1381 t459 = t31*t42
1382 t460 = t459*gamma
1383 t463 = t40*t51
1384 t464 = t463*gamma
1385 t467 = t49*t61
1386 t468 = t467*gamma
1387 t471 = t59*t185
1388 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 /0.729e3_dp*t57*t472
1392 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 /0.729e3_dp*t79*t472
1396 t490 = t66*t193*t488
1397 t493 = 0.4e1_dp/0.9e1_dp*t3*t456*t361
1398 t494 = t477*t103
1399 t496 = t102*r*t488
1400 t497 = t225*t496
1401 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 t477 + t490)*t95*t235
1405 t518 = t231*r*t516
1406 t520 = t84*t518/0.3e1_dp
1407 t522 = (-t455 + t477 - t490 + t493 - t494 + t497 + t520)*t106
1408 t525 = t2*t6
1409 t526 = t397*t525
1410 t527 = t134*t270
1411 t530 = t475*t83
1412 t532 = t274*t488
1413 t545 = (t455 - t477 + t490 + t493 - t494 + t497 + t520)*t117
1414 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 /0.3e1_dp
1422 e_tau = e_tau + (rho*t563*t102/0.8e1_dp - t124*t231*t516/0.24e2_dp)*sx
1423 t572 = t33*t141*t109
1424 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 *t57*t471
1428 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 *t79*t471
1432 t600 = t66*t193*t598
1433 t605 = t3*t450*t141*t109*t103/0.9e1_dp
1434 t606 = t587*t103
1435 t608 = t102*r*t598
1436 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 - t587 + t600)*t95*t235
1440 t630 = t231*r*t628
1441 t632 = t84*t630/0.3e1_dp
1442 t634 = (t574 + t587 - t600 - t605 - t606 + t609 + t632)*t106
1443 t639 = t585*t83
1444 t641 = t274*t598
1445 t645 = t525*t33
1446 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 *t114 - t301*t608 - t120*t630/0.3e1_dp
1453 e_laplace_rho = e_laplace_rho + (rho*t672*t102/0.8e1_dp - t124*t231*t628/0.24e2_dp)*sx
1454 END IF
1455
1456 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 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 t1 = pi**(0.1e1_dp/0.3e1_dp)
1503 t2 = t1**2
1504 t3 = br_a1*t2
1505 t4 = rho**(0.1e1_dp/0.3e1_dp)
1506 t5 = t4**2
1507 t6 = t5*rho
1508 t9 = ndrho**2
1509 t10 = 0.1e1_dp/rho
1510 t16 = laplace_rho/0.6e1_dp - gamma*(real(2*tau, kind=dp) - t9*t10/0.4e1_dp)/0.3e1_dp
1511 t17 = 0.1e1_dp/t16
1512 t18 = t6*t17
1513 t21 = 0.2e1_dp/0.3e1_dp*t3*t18 + br_a2
1514 t22 = atan(t21)
1515 t23 = -t22 + br_a3
1516 t24 = br_c1*t2
1517 t27 = t1*pi
1518 t28 = br_c2*t27
1519 t29 = rho**2
1520 t30 = t29*rho
1521 t31 = t4*t30
1522 t32 = t16**2
1523 t33 = 0.1e1_dp/t32
1524 t34 = t31*t33
1525 t37 = pi**2
1526 t38 = br_c3*t37
1527 t39 = t29**2
1528 t40 = t39*rho
1529 t42 = 0.1e1_dp/t32/t16
1530 t43 = t40*t42
1531 t46 = t2*t37
1532 t47 = br_c4*t46
1533 t48 = t39*t29
1534 t49 = t5*t48
1535 t50 = t32**2
1536 t51 = 0.1e1_dp/t50
1537 t52 = t49*t51
1538 t56 = t1*t37*pi
1539 t57 = br_c5*t56
1540 t58 = t39**2
1541 t59 = t4*t58
1542 t61 = 0.1e1_dp/t50/t16
1543 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 /0.243e3_dp*t57*t62
1547 t66 = t23*t65
1548 t67 = br_b1*t2
1549 t70 = br_b2*t27
1550 t73 = br_b3*t37
1551 t76 = br_b4*t46
1552 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 /0.243e3_dp*t79*t62
1556 t83 = 0.1e1_dp/t82
1557 t84 = t66*t83
1558 t85 = 8**(0.1e1_dp/0.3e1_dp)
1559 t86 = t23**2
1560 t87 = t86*t23
1561 t88 = t65**2
1562 t89 = t88*t65
1563 t90 = t87*t89
1564 t91 = t82**2
1565 t93 = 0.1e1_dp/t91/t82
1566 t94 = t90*t93
1567 t95 = exp(-t84)
1568 t96 = 0.1e1_dp/0.3141592654e1_dp
1569 t97 = t95*t96
1570 t99 = t94*t97*t10
1571 t100 = t99**(0.1e1_dp/0.3e1_dp)
1572 t101 = 0.1e1_dp/t100
1573 t102 = real(t85, kind=dp)*t101
1574 t103 = t102*r
1575 t104 = t84*t103
1576 t106 = exp(0.2e1_dp*t104)
1577 t108 = t106*r
1578 t109 = t108*t23
1579 t110 = t65*t83
1580 t111 = t110*t102
1581 t113 = t106*t23
1582 t115 = t84 + t104
1583 t116 = exp(t115)
1584 t118 = 0.2e1_dp*t106 - t109*t111 + t113*t110 + 0.2e1_dp + t84 + t104 &
1585 - 0.4e1_dp*t116
1586 t119 = rho*t118
1587 t121 = exp(-t115)
1588 t123 = t121*real(t85, kind=dp)*t101
1589 e_0 = e_0 + (t119*t123/0.8e1_dp)*sx
1590 IF (grad_deriv >= 1) THEN
1591 t128 = t5*t17
1592 t131 = 0.1e1_dp/t4
1593 t133 = t33*gamma
1594 t134 = t133*t9
1595 t137 = 0.10e2_dp/0.9e1_dp*t3*t128 + t3*t131*t134/0.18e2_dp
1596 t138 = t21**2
1597 t140 = 0.1e1_dp/(0.1e1_dp + t138)
1598 t141 = t137*t140
1599 t143 = t83*real(t85, kind=dp)
1600 t146 = t141*t65*t143*t101*r
1601 t153 = t4*t29
1602 t154 = t153*t33
1603 t157 = t4*rho
1604 t159 = t42*gamma
1605 t160 = t159*t9
1606 t163 = t39*t42
1607 t167 = t51*gamma
1608 t168 = t167*t9
1609 t171 = t5*t40
1610 t172 = t171*t51
1611 t175 = t5*t39
1612 t177 = t61*gamma
1613 t178 = t177*t9
1614 t182 = t4*t39*t30
1615 t183 = t182*t61
1616 t186 = t4*t48
1617 t189 = 0.1e1_dp/t50/t32
1618 t190 = t189*gamma
1619 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 t186*t191
1626 t196 = t23*t194*t83
1627 t197 = t196*t103
1628 t199 = 0.1e1_dp/t91
1629 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 t186*t191
1636 t229 = t200*t102*r*t226
1637 t232 = 0.1e1_dp/t100/t99
1638 t233 = real(t85, kind=dp)*t232
1639 t234 = t86*t89
1640 t235 = t93*t95
1641 t237 = t96*t10
1642 t242 = t87*t88*t93
1643 t247 = t91**2
1644 t249 = t90/t247
1645 t254 = t141*t110
1646 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 t256)*t95*t237 - t94*t97/t29
1650 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 *t267)*t106
1653 t272 = r*t23
1654 t277 = t194*t83
1655 t280 = t108*t66
1656 t281 = t199*real(t85, kind=dp)
1657 t292 = t140*t65*t83
1658 t295 = t65*t199
1659 t298 = t267/0.3e1_dp
1660 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 - 0.4e1_dp*t299*t116
1666 t310 = t119*t121
1667 e_rho = e_rho + (t118*t121*t102/0.8e1_dp + rho*t302*t123/0.8e1_dp - t119 &
1668 *t299*t123/0.8e1_dp - t310*t233*t264/0.24e2_dp)*sx
1669 t314 = t3*t5
1670 t315 = t133*ndrho
1671 t317 = t292*t103
1672 t318 = t314*t315*t317
1673 t324 = t159*ndrho
1674 t328 = t167*ndrho
1675 t332 = t177*ndrho
1676 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 *t332 - 0.80e2_dp/0.729e3_dp*t57*t182*t336
1680 t341 = t23*t339*t83
1681 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 *t332 - 0.80e2_dp/0.729e3_dp*t79*t182*t336
1685 t362 = t200*t102*r*t359
1686 t368 = gamma*ndrho
1687 t369 = t368*t140
1688 t383 = t368*t292
1689 t385 = t3*t5*t33*t383/0.9e1_dp
1690 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 (-t385 - t341 + t387)*t95*t237
1694 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 /0.3e1_dp*t395)*t106
1697 t402 = t108*br_a1
1698 t404 = t2*t5*t33
1699 t409 = t339*t83
1700 t420 = t106*br_a1
1701 t427 = t318/0.9e1_dp
1702 t428 = t395/0.3e1_dp
1703 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 + 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 *t233*t392/0.24e2_dp)*sx
1711 t443 = t6*t33
1712 t444 = t443*gamma
1713 t446 = t3*t444*t317
1714 t450 = t31*t42
1715 t451 = t450*gamma
1716 t454 = t40*t51
1717 t455 = t454*gamma
1718 t458 = t49*t61
1719 t459 = t458*gamma
1720 t462 = t59*t189
1721 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 /0.729e3_dp*t57*t463
1725 t468 = t23*t466*t83
1726 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 /0.729e3_dp*t79*t463
1730 t484 = t200*t102*r*t481
1731 t487 = t234*t235*t96
1732 t501 = gamma*t140
1733 t504 = 0.4e1_dp/0.9e1_dp*t3*t443*t501*t110
1734 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 t468 + t506)*t95*t237
1738 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 /0.3e1_dp*t514)*t106
1741 t521 = t2*t6
1742 t525 = t143*t101
1743 t529 = t466*t83
1744 t540 = t420*t521
1745 t547 = 0.4e1_dp/0.9e1_dp*t446
1746 t548 = t514/0.3e1_dp
1747 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 *t116
1754 e_tau = e_tau + (rho*t552*t123/0.8e1_dp - t119*t549*t123/0.8e1_dp - t310 &
1755 *t233*t511/0.24e2_dp)*sx
1756 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 *t57*t462
1760 t580 = t23*t578*t83
1761 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 *t79*t462
1765 t596 = t200*t102*r*t593
1766 t612 = t3*t6
1767 t613 = t33*t140
1768 t614 = t613*t110
1769 t616 = t612*t614/0.9e1_dp
1770 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 - t580 + t618)*t95*t237
1774 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 /0.3e1_dp*t626)*t106
1777 t638 = t578*t83
1778 t654 = t566/0.9e1_dp
1779 t655 = t626/0.3e1_dp
1780 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 + 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 *t233*t623/0.24e2_dp)*sx
1788 END IF
1789
1790 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 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 t1 = 8**(0.1e1_dp/0.3e1_dp)
1824 t2 = t1**2
1825 t3 = 0.1e1_dp/br_bb
1826 t4 = pi**(0.1e1_dp/0.3e1_dp)
1827 t5 = t4**2
1828 t6 = 0.1e1_dp/t5
1829 t8 = rho**(0.1e1_dp/0.3e1_dp)
1830 t9 = t8**2
1831 t10 = t9*rho
1832 t11 = 0.1e1_dp/t10
1833 t14 = ndrho**2
1834 t15 = 0.1e1_dp/rho
1835 t21 = laplace_rho/0.6e1_dp - gamma*(real(2*tau, kind=dp) - t14*t15/0.4e1_dp)/0.3e1_dp
1836 t25 = br_bb**2
1837 t26 = t4*pi
1838 t28 = rho**2
1839 t31 = t21**2
1840 t33 = t8*t28*rho/t31
1841 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 *t37*t3*t6*t11*t21)
1844 t45 = t44 + 0.2e1_dp
1845 t46 = t45**2
1846 t50 = t10/t21
1847 t56 = pi**2
1848 t58 = t28**2
1849 t62 = t58*rho/t31/t21
1850 t65 = t5*t56
1851 t69 = t31**2
1852 t71 = t9*t58*t28/t69
1853 t75 = t4*t56*pi
1854 t77 = t58**2
1855 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 *br_d4*t65*t71 + 0.32e2_dp/0.243e3_dp*br_d5*t75*t81
1859 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 *br_e4*t65*t71 + 0.32e2_dp/0.243e3_dp*br_e5*t75*t81
1863 t104 = t103**2
1864 t111 = exp(-t45*t84/t103)
1865 t116 = (t46*t45*t85*t84/t104/t103*t111/0.3141592654e1_dp &
1866 *t15)**(0.1e1_dp/0.3e1_dp)
1867 brval = real(t2, kind=dp)*t116/0.8e1_dp
1868
1869 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 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 sx, r, gamma, grad_deriv)
1877 END IF
1878
1879 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 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 t1 = 0.1e1_dp/br_bb
1927 t2 = pi**(0.1e1_dp/0.3e1_dp)
1928 t3 = t2**2
1929 t4 = 0.1e1_dp/t3
1930 t5 = t1*t4
1931 t6 = rho**(0.1e1_dp/0.3e1_dp)
1932 t7 = t6**2
1933 t8 = t7*rho
1934 t9 = 0.1e1_dp/t8
1935 t12 = ndrho**2
1936 t13 = 0.1e1_dp/rho
1937 t19 = laplace_rho/0.6e1_dp - gamma*(real(2*tau, kind=dp) - t12*t13/0.4e1_dp)/0.3e1_dp
1938 t20 = t9*t19
1939 t23 = br_bb**2
1940 t24 = t2*pi
1941 t25 = t23*t24
1942 t26 = rho**2
1943 t27 = t26*rho
1944 t28 = t6*t27
1945 t29 = t19**2
1946 t30 = 0.1e1_dp/t29
1947 t31 = t28*t30
1948 t35 = sqrt(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t31)
1949 t36 = t35*t1
1950 t37 = t4*t9
1951 t41 = 0.1500000000000000e1_dp*t5*t20 + 0.3e1_dp/0.2e1_dp*t36*t37* &
1952 t19
1953 t42 = log(t41)
1954 t43 = t42 + 0.2e1_dp
1955 t44 = br_d1*t3
1956 t45 = 0.1e1_dp/t19
1957 t46 = t8*t45
1958 t49 = br_d2*t24
1959 t52 = pi**2
1960 t53 = br_d3*t52
1961 t54 = t26**2
1962 t55 = t54*rho
1963 t57 = 0.1e1_dp/t29/t19
1964 t58 = t55*t57
1965 t61 = t3*t52
1966 t62 = br_d4*t61
1967 t63 = t54*t26
1968 t64 = t7*t63
1969 t65 = t29**2
1970 t66 = 0.1e1_dp/t65
1971 t67 = t64*t66
1972 t71 = t2*t52*pi
1973 t72 = br_d5*t71
1974 t73 = t54**2
1975 t74 = t6*t73
1976 t76 = 0.1e1_dp/t65/t19
1977 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 /0.243e3_dp*t72*t77
1981 t81 = t43*t80
1982 t82 = br_e1*t3
1983 t85 = br_e2*t24
1984 t88 = br_e3*t52
1985 t91 = br_e4*t61
1986 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 /0.243e3_dp*t94*t77
1990 t98 = 0.1e1_dp/t97
1991 t99 = t81*t98
1992 t100 = 8**(0.1e1_dp/0.3e1_dp)
1993 t101 = t43**2
1994 t102 = t101*t43
1995 t103 = t80**2
1996 t104 = t103*t80
1997 t105 = t102*t104
1998 t106 = t97**2
1999 t108 = 0.1e1_dp/t106/t97
2000 t109 = t105*t108
2001 t110 = exp(-t99)
2002 t111 = 0.1e1_dp/0.3141592654e1_dp
2003 t112 = t110*t111
2004 t114 = t109*t112*t13
2005 t115 = t114**(0.1e1_dp/0.3e1_dp)
2006 t116 = 0.1e1_dp/t115
2007 t117 = real(t100, kind=dp)*t116
2008 t118 = t117*r
2009 t119 = t99*t118
2010 t121 = exp(t99 - t119)
2011 t123 = t121*t43
2012 t124 = t80*t98
2013 t125 = t123*t124
2014 t129 = t98*real(t100, kind=dp)*t116*r
2015 t132 = exp(-t99 - t119)
2016 t134 = t132*t43
2017 t135 = t134*t124
2018 t138 = -0.2e1_dp*t121 + t125 - t123*t80*t129 + 0.2e1_dp*t132 + t135 &
2019 + t134*t80*t129
2020 t139 = rho*t138
2021 e_0 = e_0 + (t139*t117/0.8e1_dp)*sx
2022 IF (grad_deriv >= 1) THEN
2023 t145 = 0.1e1_dp/t7/t26
2024 t152 = 0.1e1_dp/t7/t27*gamma*t12
2025 t155 = 0.1e1_dp/t35
2026 t158 = t6*t26
2027 t159 = t158*t30
2028 t162 = t6*rho
2029 t164 = t57*gamma
2030 t165 = t164*t12
2031 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 *t4*t145*t19 - t176*t152/0.8e1_dp
2036 t180 = 0.1e1_dp/t41
2037 t181 = t179*t180
2038 t182 = t181*t124
2039 t183 = t7*t45
2040 t186 = 0.1e1_dp/t6
2041 t188 = t30*gamma
2042 t189 = t188*t12
2043 t197 = t54*t57
2044 t201 = t66*gamma
2045 t202 = t201*t12
2046 t205 = t7*t55
2047 t206 = t205*t66
2048 t209 = t7*t54
2049 t211 = t76*gamma
2050 t212 = t211*t12
2051 t216 = t6*t54*t27
2052 t217 = t216*t76
2053 t220 = t6*t63
2054 t223 = 0.1e1_dp/t65/t29
2055 t224 = t223*gamma
2056 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 t220*t225
2063 t230 = t43*t228*t98
2064 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 t220*t225
2071 t259 = t81*t231*t257
2072 t261 = t181*t80*t129
2073 t262 = t230*t118
2074 t263 = t81*t231
2075 t265 = t117*r*t257
2076 t266 = t263*t265
2077 t269 = real(t100, kind=dp)/t115/t114
2078 t272 = t101*t104*t108*t110
2079 t273 = t111*t13
2080 t278 = t102*t103*t108
2081 t283 = t106**2
2082 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 *t110*t273 - t109*t112/t26
2086 t299 = t269*r*t297
2087 t301 = t99*t299/0.3e1_dp
2088 t303 = (t182 + t230 - t259 - t261 - t262 + t266 + t301)*t121
2089 t306 = t121*t179
2090 t307 = t180*t80
2091 t308 = t307*t98
2092 t310 = t228*t98
2093 t312 = t80*t231
2094 t313 = t312*t257
2095 t321 = t123*t312
2096 t326 = (-t182 - t230 + t259 - t261 - t262 + t266 + t301)*t132
2097 t329 = t132*t179
2098 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 /0.3e1_dp
2105 e_rho = e_rho + (t138*real(t100, kind=dp)*t116/0.8e1_dp + rho*t343*t117/0.8e1_dp &
2106 - t139*t269*t297/0.24e2_dp)*sx
2107 t351 = t145*gamma*ndrho
2108 t354 = t155*br_bb
2109 t355 = t354*t3
2110 t363 = 0.2500000000000000e0_dp*t5*t351 - t355*t7*t30*gamma*ndrho &
2111 /0.9e1_dp + t176*t351/0.4e1_dp
2112 t364 = t363*t180
2113 t365 = t364*t124
2114 t367 = t188*ndrho
2115 t371 = t164*ndrho
2116 t375 = t201*ndrho
2117 t379 = t211*ndrho
2118 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 *t379 - 0.80e2_dp/0.729e3_dp*t72*t216*t383
2122 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 *t379 - 0.80e2_dp/0.729e3_dp*t94*t216*t383
2126 t406 = t81*t231*t404
2127 t408 = t364*t80*t129
2128 t409 = t388*t118
2129 t411 = t117*r*t404
2130 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 *t110*t273
2134 t430 = t269*r*t428
2135 t432 = t99*t430/0.3e1_dp
2136 t434 = (t365 + t388 - t406 - t408 - t409 + t412 + t432)*t121
2137 t437 = t121*t363
2138 t439 = t386*t98
2139 t441 = t312*t404
2140 t453 = (-t365 - t388 + t406 - t408 - t409 + t412 + t432)*t132
2141 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 /0.3e1_dp
2148 e_ndrho = e_ndrho + (rho*t469*t117/0.8e1_dp - t139*t269*t428/0.24e2_dp)*sx
2149 t479 = t8*t30
2150 t480 = t479*gamma
2151 t485 = -0.1000000000e1_dp*t5*t9*gamma + 0.4e1_dp/0.9e1_dp*t355*t480 &
2152 - t36*t37*gamma
2153 t486 = t485*t180
2154 t487 = t486*t124
2155 t490 = t28*t57
2156 t491 = t490*gamma
2157 t494 = t55*t66
2158 t495 = t494*gamma
2159 t498 = t64*t76
2160 t499 = t498*gamma
2161 t502 = t74*t223
2162 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 /0.729e3_dp*t72*t503
2166 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 /0.729e3_dp*t94*t503
2170 t521 = t81*t231*t519
2171 t523 = t486*t80*t129
2172 t524 = t508*t118
2173 t526 = t117*r*t519
2174 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 *t110*t273
2178 t545 = t269*r*t543
2179 t547 = t99*t545/0.3e1_dp
2180 t549 = (t487 + t508 - t521 - t523 - t524 + t527 + t547)*t121
2181 t552 = t121*t485
2182 t554 = t506*t98
2183 t556 = t312*t519
2184 t568 = (-t487 - t508 + t521 - t523 - t524 + t527 + t547)*t132
2185 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 /0.3e1_dp
2192 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 + t36*t37/0.4e1_dp
2195 t600 = t599*t180
2196 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 *t72*t502
2200 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 *t94*t502
2204 t627 = t81*t231*t625
2205 t629 = t600*t80*t129
2206 t630 = t614*t118
2207 t632 = t117*r*t625
2208 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 *t110*t273
2212 t651 = t269*r*t649
2213 t653 = t99*t651/0.3e1_dp
2214 t655 = (t601 + t614 - t627 - t629 - t630 + t633 + t653)*t121
2215 t658 = t121*t599
2216 t660 = t612*t98
2217 t662 = t312*t625
2218 t674 = (-t601 - t614 + t627 - t629 - t630 + t633 + t653)*t132
2219 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 /0.3e1_dp
2226 e_laplace_rho = e_laplace_rho + (rho*t690*t117/0.8e1_dp - t139*t269*t649/0.24e2_dp)*sx
2227 END IF
2228
2229 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 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 t1 = 0.1e1_dp/br_bb
2276 t2 = pi**(0.1e1_dp/0.3e1_dp)
2277 t3 = t2**2
2278 t4 = 0.1e1_dp/t3
2279 t5 = t1*t4
2280 t6 = rho**(0.1e1_dp/0.3e1_dp)
2281 t7 = t6**2
2282 t8 = t7*rho
2283 t9 = 0.1e1_dp/t8
2284 t12 = ndrho**2
2285 t13 = 0.1e1_dp/rho
2286 t19 = laplace_rho/0.6e1_dp - gamma*(real(2*tau, kind=dp) - t12*t13/0.4e1_dp)/0.3e1_dp
2287 t20 = t9*t19
2288 t23 = br_bb**2
2289 t24 = t2*pi
2290 t25 = t23*t24
2291 t26 = rho**2
2292 t27 = t26*rho
2293 t28 = t6*t27
2294 t29 = t19**2
2295 t30 = 0.1e1_dp/t29
2296 t31 = t28*t30
2297 t35 = sqrt(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t31)
2298 t36 = t35*t1
2299 t37 = t4*t9
2300 t41 = 0.1500000000000000e1_dp*t5*t20 + 0.3e1_dp/0.2e1_dp*t36*t37* &
2301 t19
2302 t42 = log(t41)
2303 t43 = t42 + 0.2e1_dp
2304 t44 = br_d1*t3
2305 t45 = 0.1e1_dp/t19
2306 t46 = t8*t45
2307 t49 = br_d2*t24
2308 t52 = pi**2
2309 t53 = br_d3*t52
2310 t54 = t26**2
2311 t55 = t54*rho
2312 t57 = 0.1e1_dp/t29/t19
2313 t58 = t55*t57
2314 t61 = t3*t52
2315 t62 = br_d4*t61
2316 t63 = t54*t26
2317 t64 = t7*t63
2318 t65 = t29**2
2319 t66 = 0.1e1_dp/t65
2320 t67 = t64*t66
2321 t71 = t2*t52*pi
2322 t72 = br_d5*t71
2323 t73 = t54**2
2324 t74 = t6*t73
2325 t76 = 0.1e1_dp/t65/t19
2326 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 /0.243e3_dp*t72*t77
2330 t81 = t43*t80
2331 t82 = br_e1*t3
2332 t85 = br_e2*t24
2333 t88 = br_e3*t52
2334 t91 = br_e4*t61
2335 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 /0.243e3_dp*t94*t77
2339 t98 = 0.1e1_dp/t97
2340 t99 = t81*t98
2341 t100 = 8**(0.1e1_dp/0.3e1_dp)
2342 t101 = t43**2
2343 t102 = t101*t43
2344 t103 = t80**2
2345 t104 = t103*t80
2346 t105 = t102*t104
2347 t106 = t97**2
2348 t108 = 0.1e1_dp/t106/t97
2349 t109 = t105*t108
2350 t110 = exp(-t99)
2351 t111 = 0.1e1_dp/0.3141592654e1_dp
2352 t112 = t110*t111
2353 t114 = t109*t112*t13
2354 t115 = t114**(0.1e1_dp/0.3e1_dp)
2355 t116 = 0.1e1_dp/t115
2356 t117 = real(t100, kind=dp)*t116
2357 t118 = t117*r
2358 t119 = t99*t118
2359 t121 = exp(0.2e1_dp*t119)
2360 t123 = t121*r
2361 t124 = t123*t43
2362 t125 = t80*t98
2363 t126 = t125*t117
2364 t128 = t121*t43
2365 t130 = t99 + t119
2366 t131 = exp(t130)
2367 t133 = 0.2e1_dp*t121 - t124*t126 + t128*t125 + 0.2e1_dp + t99 + t119 &
2368 - 0.4e1_dp*t131
2369 t134 = rho*t133
2370 t136 = exp(-t130)
2371 t138 = t136*real(t100, kind=dp)*t116
2372 e_0 = e_0 + (t134*t138/0.8e1_dp)*sx
2373 IF (grad_deriv >= 1) THEN
2374 t144 = 0.1e1_dp/t7/t26
2375 t151 = 0.1e1_dp/t7/t27*gamma*t12
2376 t154 = 0.1e1_dp/t35
2377 t157 = t6*t26
2378 t158 = t157*t30
2379 t161 = t6*rho
2380 t163 = t57*gamma
2381 t164 = t163*t12
2382 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 *t4*t144*t19 - t175*t151/0.8e1_dp
2387 t179 = 0.1e1_dp/t41
2388 t180 = t178*t179
2389 t182 = t98*real(t100, kind=dp)
2390 t184 = t182*t116*r
2391 t185 = t180*t80*t184
2392 t187 = t7*t45
2393 t190 = 0.1e1_dp/t6
2394 t192 = t30*gamma
2395 t193 = t192*t12
2396 t201 = t54*t57
2397 t205 = t66*gamma
2398 t206 = t205*t12
2399 t209 = t7*t55
2400 t210 = t209*t66
2401 t213 = t7*t54
2402 t215 = t76*gamma
2403 t216 = t215*t12
2404 t220 = t6*t54*t27
2405 t221 = t220*t76
2406 t224 = t6*t63
2407 t227 = 0.1e1_dp/t65/t29
2408 t228 = t227*gamma
2409 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 t224*t229
2416 t234 = t43*t232*t98
2417 t235 = t234*t118
2418 t237 = 0.1e1_dp/t106
2419 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 t224*t229
2426 t267 = t238*t117*r*t264
2427 t270 = 0.1e1_dp/t115/t114
2428 t271 = real(t100, kind=dp)*t270
2429 t274 = t101*t104*t108*t110
2430 t275 = t111*t13
2431 t280 = t102*t103*t108
2432 t285 = t106**2
2433 t287 = t105/t285
2434 t292 = t180*t125
2435 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 *t110*t275 - t109*t112/t26
2439 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 *t305)*t121
2442 t310 = r*t43
2443 t315 = t232*t98
2444 t318 = t123*t81
2445 t319 = t237*real(t100, kind=dp)
2446 t330 = t179*t80*t98
2447 t333 = t80*t237
2448 t336 = t305/0.3e1_dp
2449 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 - 0.4e1_dp*t337*t131
2455 t348 = t134*t136
2456 e_rho = e_rho + (t133*t136*t117/0.8e1_dp + rho*t340*t138/0.8e1_dp - t134 &
2457 *t337*t138/0.8e1_dp - t348*t271*t302/0.24e2_dp)*sx
2458 t353 = t144*gamma*ndrho
2459 t356 = t154*br_bb
2460 t357 = t356*t3
2461 t365 = 0.2500000000000000e0_dp*t5*t353 - t357*t7*t30*gamma*ndrho &
2462 /0.9e1_dp + t175*t353/0.4e1_dp
2463 t366 = t365*t179
2464 t368 = t366*t80*t184
2465 t371 = t192*ndrho
2466 t375 = t163*ndrho
2467 t379 = t205*ndrho
2468 t383 = t215*ndrho
2469 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 *t383 - 0.80e2_dp/0.729e3_dp*t72*t220*t387
2473 t392 = t43*t390*t98
2474 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 *t383 - 0.80e2_dp/0.729e3_dp*t94*t220*t387
2478 t413 = t238*t117*r*t410
2479 t426 = t366*t125
2480 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 *t110*t275
2484 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 *t436)*t121
2487 t445 = t390*t98
2488 t461 = t436/0.3e1_dp
2489 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 - 0.4e1_dp*t462*t131
2495 e_ndrho = e_ndrho + (rho*t465*t138/0.8e1_dp - t134*t462*t138/0.8e1_dp - t348 &
2496 *t271*t433/0.24e2_dp)*sx
2497 t479 = t8*t30
2498 t480 = t479*gamma
2499 t485 = -0.1000000000e1_dp*t5*t9*gamma + 0.4e1_dp/0.9e1_dp*t357*t480 &
2500 - t36*t37*gamma
2501 t486 = t485*t179
2502 t488 = t486*t80*t184
2503 t492 = t28*t57
2504 t493 = t492*gamma
2505 t496 = t55*t66
2506 t497 = t496*gamma
2507 t500 = t64*t76
2508 t501 = t500*gamma
2509 t504 = t74*t227
2510 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 /0.729e3_dp*t72*t505
2514 t510 = t43*t508*t98
2515 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 /0.729e3_dp*t94*t505
2519 t526 = t238*t117*r*t523
2520 t539 = t486*t125
2521 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 *t110*t275
2525 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 *t549)*t121
2528 t558 = t508*t98
2529 t574 = t549/0.3e1_dp
2530 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 - 0.4e1_dp*t575*t131
2536 e_tau = e_tau + (rho*t578*t138/0.8e1_dp - t134*t575*t138/0.8e1_dp - t348 &
2537 *t271*t546/0.24e2_dp)*sx
2538 t597 = 0.2500000000000000e0_dp*t5*t9 - t356*t3*t8*t30/0.9e1_dp &
2539 + t36*t37/0.4e1_dp
2540 t598 = t597*t179
2541 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 *t72*t504
2545 t614 = t43*t612*t98
2546 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 *t94*t504
2550 t630 = t238*t117*r*t627
2551 t643 = t598*t125
2552 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 *t110*t275
2556 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 *t653)*t121
2559 t662 = t612*t98
2560 t678 = t653/0.3e1_dp
2561 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 - 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 *t271*t650/0.24e2_dp)*sx
2569 END IF
2570
2571 END SUBROUTINE x_br_lsd_y_gt_0_cutoff_r_lte_b
2572
2573END MODULE xc_xbecke_roussel
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::proynov2007
integer, save, public proynov2007
Definition bibliography.F:36
bibliography::beckeroussel1989
integer, save, public beckeroussel1989
Definition bibliography.F:36
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition xc_derivative_desc.F:26
xc_derivative_desc::deriv_norm_drho
integer, parameter, public deriv_norm_drho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_laplace_rhob
integer, parameter, public deriv_laplace_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau
integer, parameter, public deriv_tau
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau_b
integer, parameter, public deriv_tau_b
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau_a
integer, parameter, public deriv_tau_a
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_laplace_rhoa
integer, parameter, public deriv_laplace_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_laplace_rho
integer, parameter, public deriv_laplace_rho
Definition xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:65
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition xc_rho_set_types.F:281
xc_xbecke_roussel
Calculates the exchange energy based on the Becke-Roussel exchange hole. Takes advantage of an analyt...
Definition xc_xbecke_roussel.F:18
xc_xbecke_roussel::x_br_lsd_y_gt_0
subroutine, public x_br_lsd_y_gt_0(rho, ndrho, tau, laplace_rho, e_0, e_rho, e_ndrho, e_tau, e_laplace_rho, sx, gamma, grad_deriv)
Full range evaluation for y > 0.
Definition xc_xbecke_roussel.F:804
xc_xbecke_roussel::x_br_lsd_y_gt_0_cutoff
subroutine, public x_br_lsd_y_gt_0_cutoff(rho, ndrho, tau, laplace_rho, e_0, e_rho, e_ndrho, e_tau, e_laplace_rho, sx, r, gamma, grad_deriv)
Truncated long range part for y > 0.
Definition xc_xbecke_roussel.F:1814
xc_xbecke_roussel::xbecke_roussel_lsd_eval
subroutine, public xbecke_roussel_lsd_eval(rho_set, deriv_set, grad_deriv, br_params)
evaluates the Becke Roussel exchange functional for lda
Definition xc_xbecke_roussel.F:353
xc_xbecke_roussel::x_br_lsd_y_lte_0_cutoff
subroutine, public x_br_lsd_y_lte_0_cutoff(rho, ndrho, tau, laplace_rho, e_0, e_rho, e_ndrho, e_tau, e_laplace_rho, sx, r, gamma, grad_deriv)
Truncated long range part for y <= 0.
Definition xc_xbecke_roussel.F:1060
xc_xbecke_roussel::xbecke_roussel_lda_info
subroutine, public xbecke_roussel_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_xbecke_roussel.F:91
xc_xbecke_roussel::x_br_lsd_y_lte_0
subroutine, public x_br_lsd_y_lte_0(rho, ndrho, tau, laplace_rho, e_0, e_rho, e_ndrho, e_tau, e_laplace_rho, sx, gamma, grad_deriv)
full range evaluation for y <= 0
Definition xc_xbecke_roussel.F:562
xc_xbecke_roussel::xbecke_roussel_lda_eval
subroutine, public xbecke_roussel_lda_eval(rho_set, deriv_set, grad_deriv, br_params)
evaluates the Becke Roussel exchange functional for lda
Definition xc_xbecke_roussel.F:162
xc_xbecke_roussel::xbecke_roussel_lsd_info
subroutine, public xbecke_roussel_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_xbecke_roussel.F:127
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78