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