LCOV - code coverage report
Current view: top level - src/xc - xc_xbecke_roussel.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.9 % 1512 1510
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 14 14

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Calculates the exchange energy based on the Becke-Roussel exchange
      10              : !>        hole. Takes advantage of an analytical representation of the hole
      11              : !>        in order to avoid solving a non-linear equation by means of Newton-
      12              : !>        Raphson algorithm
      13              : !> \note
      14              : !> \par History
      15              : !>      12.2008 created [mguidon]
      16              : !> \author mguidon
      17              : ! **************************************************************************************************
      18              : MODULE xc_xbecke_roussel
      19              :    USE bibliography,                    ONLY: BeckeRoussel1989,&
      20              :                                               Proynov2007,&
      21              :                                               cite_reference
      22              :    USE input_section_types,             ONLY: section_vals_type,&
      23              :                                               section_vals_val_get
      24              :    USE kinds,                           ONLY: dp
      25              :    USE mathconstants,                   ONLY: pi
      26              :    USE xc_derivative_desc,              ONLY: &
      27              :         deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, deriv_norm_drho, &
      28              :         deriv_norm_drhoa, deriv_norm_drhob, deriv_rho, deriv_rhoa, deriv_rhob, deriv_tau, &
      29              :         deriv_tau_a, deriv_tau_b
      30              :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      31              :                                               xc_dset_get_derivative
      32              :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      33              :                                               xc_derivative_type
      34              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      35              :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      36              :                                               xc_rho_set_type
      37              : #include "../base/base_uses.f90"
      38              : 
      39              :    IMPLICIT NONE
      40              :    PRIVATE
      41              : 
      42              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      43              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbecke_roussel'
      44              : 
      45              :    REAL(dp), PARAMETER, PRIVATE  :: br_a1 = 1.5255251812009530_dp, &
      46              :                                     br_a2 = 0.4576575543602858_dp, &
      47              :                                     br_a3 = 0.4292036732051034_dp, &
      48              :                                     br_c0 = 0.7566445420735584_dp, &
      49              :                                     br_c1 = -2.6363977871370960_dp, &
      50              :                                     br_c2 = 5.4745159964232880_dp, &
      51              :                                     br_c3 = -12.657308127108290_dp, &
      52              :                                     br_c4 = 4.1250584725121360_dp, &
      53              :                                     br_c5 = -30.425133957163840_dp, &
      54              :                                     br_b0 = 0.4771976183772063_dp, &
      55              :                                     br_b1 = -1.7799813494556270_dp, &
      56              :                                     br_b2 = 3.8433841862302150_dp, &
      57              :                                     br_b3 = -9.5912050880518490_dp, &
      58              :                                     br_b4 = 2.1730180285916720_dp, &
      59              :                                     br_b5 = -30.425133851603660_dp, &
      60              :                                     br_d0 = 0.00004435009886795587_dp, &
      61              :                                     br_d1 = 0.58128653604457910_dp, &
      62              :                                     br_d2 = 66.742764515940610_dp, &
      63              :                                     br_d3 = 434.26780897229770_dp, &
      64              :                                     br_d4 = 824.7765766052239000_dp, &
      65              :                                     br_d5 = 1657.9652731582120_dp, &
      66              :                                     br_e0 = 0.00003347285060926091_dp, &
      67              :                                     br_e1 = 0.47917931023971350_dp, &
      68              :                                     br_e2 = 62.392268338574240_dp, &
      69              :                                     br_e3 = 463.14816427938120_dp, &
      70              :                                     br_e4 = 785.2360350104029000_dp, &
      71              :                                     br_e5 = 1657.962968223273000000_dp, &
      72              :                                     br_BB = 2.085749716493756_dp
      73              : 
      74              :    PUBLIC :: xbecke_roussel_lda_info, xbecke_roussel_lsd_info, xbecke_roussel_lda_eval, &
      75              :              xbecke_roussel_lsd_eval, x_br_lsd_y_lte_0, x_br_lsd_y_gt_0, x_br_lsd_y_lte_0_cutoff, &
      76              :              x_br_lsd_y_gt_0_cutoff
      77              : CONTAINS
      78              : 
      79              : ! **************************************************************************************************
      80              : !> \brief return various information on the functional
      81              : !> \param reference string with the reference of the actual functional
      82              : !> \param shortform string with the shortform of the functional name
      83              : !> \param needs the components needed by this functional are set to
      84              : !>        true (does not set the unneeded components to false)
      85              : !> \param max_deriv controls the number of derivatives
      86              : !> \par History
      87              : !>        12.2008 created [mguidon]
      88              : !> \author mguidon
      89              : ! **************************************************************************************************
      90           99 :    SUBROUTINE xbecke_roussel_lda_info(reference, shortform, needs, max_deriv)
      91              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      92              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      93              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      94              : 
      95           99 :       CALL cite_reference(BeckeRoussel1989)
      96           99 :       CALL cite_reference(Proynov2007)
      97           99 :       IF (PRESENT(reference)) THEN
      98              :          reference = "A.D. Becke, M.R. Roussel, "// &
      99            3 :                      "Phys. Rev. A, vol. 39, n 8, pp. 3761-3767, (1989) {spin unpolarized}"
     100              :       END IF
     101           99 :       IF (PRESENT(shortform)) THEN
     102            3 :          shortform = "Becke-Roussel exchange hole (spin unpolarized)"
     103              :       END IF
     104           99 :       IF (PRESENT(needs)) THEN
     105           96 :          needs%rho = .TRUE.
     106           96 :          needs%norm_drho = .TRUE.
     107           96 :          needs%tau = .TRUE.
     108           96 :          needs%laplace_rho = .TRUE.
     109              :       END IF
     110              : 
     111           99 :       IF (PRESENT(max_deriv)) max_deriv = 1
     112              : 
     113           99 :    END SUBROUTINE xbecke_roussel_lda_info
     114              : 
     115              : ! **************************************************************************************************
     116              : !> \brief return various information on the functional
     117              : !> \param reference string with the reference of the actual functional
     118              : !> \param shortform string with the shortform of the functional name
     119              : !> \param needs the components needed by this functional are set to
     120              : !>        true (does not set the unneeded components to false)
     121              : !> \param max_deriv ...
     122              : !> \par History
     123              : !>        12.2008 created [mguidon]
     124              : !> \author mguidon
     125              : ! **************************************************************************************************
     126           92 :    SUBROUTINE xbecke_roussel_lsd_info(reference, shortform, needs, max_deriv)
     127              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
     128              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
     129              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
     130              : 
     131           92 :       CALL cite_reference(BeckeRoussel1989)
     132           92 :       CALL cite_reference(Proynov2007)
     133           92 :       IF (PRESENT(reference)) THEN
     134              :          reference = "A.D. Becke, M.R. Roussel, "// &
     135            2 :                      "Phys. Rev. A, vol. 39, n 8, pp. 3761-3767, (1989) {spin polarized}"
     136              :       END IF
     137           92 :       IF (PRESENT(shortform)) THEN
     138            2 :          shortform = "Becke-Roussel exchange hole (spin polarized)"
     139              :       END IF
     140           92 :       IF (PRESENT(needs)) THEN
     141           90 :          needs%rho_spin = .TRUE.
     142           90 :          needs%norm_drho_spin = .TRUE.
     143           90 :          needs%tau_spin = .TRUE.
     144           90 :          needs%laplace_rho_spin = .TRUE.
     145              :       END IF
     146           92 :       IF (PRESENT(max_deriv)) max_deriv = 1
     147              : 
     148           92 :    END SUBROUTINE xbecke_roussel_lsd_info
     149              : 
     150              : ! **************************************************************************************************
     151              : !> \brief evaluates the Becke Roussel exchange functional for lda
     152              : !> \param rho_set the density where you want to evaluate the functional
     153              : !> \param deriv_set place where to store the functional derivatives (they are
     154              : !>        added to the derivatives)
     155              : !> \param grad_deriv 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          420 :    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           84 :          POINTER                                         :: dummy, e_0, e_laplace_rho, e_ndrho, &
     177           84 :                                                             e_rho, e_tau, laplace_rho, norm_drho, &
     178           84 :                                                             rho, tau
     179              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     180              : 
     181           84 :       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           84 :                           rho_cutoff=epsilon_rho)
     186           84 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     187              : 
     188           84 :       dummy => rho
     189              : 
     190           84 :       e_0 => dummy
     191           84 :       e_rho => dummy
     192           84 :       e_ndrho => dummy
     193           84 :       e_tau => dummy
     194           84 :       e_laplace_rho => dummy
     195              : 
     196           84 :       IF (grad_deriv >= 0) THEN
     197              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     198           84 :                                          allocate_deriv=.TRUE.)
     199           84 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     200              :       END IF
     201           84 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     202              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     203           84 :                                          allocate_deriv=.TRUE.)
     204           84 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     205              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     206           84 :                                          allocate_deriv=.TRUE.)
     207           84 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     208              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
     209           84 :                                          allocate_deriv=.TRUE.)
     210           84 :          CALL xc_derivative_get(deriv, deriv_data=e_tau)
     211              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho], &
     212           84 :                                          allocate_deriv=.TRUE.)
     213           84 :          CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho)
     214              :       END IF
     215           84 :       IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
     216            0 :          CPABORT("derivatives bigger than 1 not implemented")
     217              :       END IF
     218              : 
     219           84 :       CALL section_vals_val_get(br_params, "scale_x", r_val=sx)
     220           84 :       CALL section_vals_val_get(br_params, "CUTOFF_RADIUS", r_val=R)
     221           84 :       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           84 : !$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           84 :       CALL timestop(handle)
     238           84 :    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           84 :    SUBROUTINE xbecke_roussel_lda_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
     268           84 :                                       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      1283968 :          my_rho = 0.5_dp*MAX(rho(ip), 0.0_dp)
     287      1283968 :          IF (my_rho > epsilon_rho) THEN
     288      1283968 :             my_ndrho = 0.5_dp*MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp)
     289      1283968 :             my_tau = 0.5_dp*MAX(EPSILON(0.0_dp)*1.e4_dp, tau(ip))
     290      1283968 :             my_laplace_rho = 0.5_dp*laplace_rho(ip)
     291              : 
     292      1283968 :             t1 = pi**(0.1e1_dp/0.3e1_dp)
     293      1283968 :             t2 = t1**2
     294      1283968 :             t3 = my_rho**(0.1e1_dp/0.3e1_dp)
     295      1283968 :             t4 = t3**2
     296      1283968 :             t5 = t4*my_rho
     297      1283968 :             t8 = my_ndrho**2
     298      1283968 :             t9 = 0.1e1_dp/my_rho
     299              :             ! *** CP2K defines tau in a different way as compared to Becke !!!
     300      1283968 :             t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
     301      1283968 :             t16 = 0.1e1_dp/t15
     302      1283968 :             yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
     303      1283968 :             IF (R == 0.0_dp) THEN
     304       890752 :                IF (yval <= 0.0_dp) THEN
     305       138801 :                   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       138801 :                                         sx, gamma, grad_deriv)
     309              :                   ! VERY UGLY HACK e_0 has to multiplied by the factor 2
     310       138801 :                   e_diff = e_0(ip) - e_old
     311       138801 :                   e_0(ip) = e_0(ip) + e_diff
     312              :                ELSE
     313       751951 :                   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       751951 :                                        sx, gamma, grad_deriv)
     317              :                   ! VERY UGLY HACK e_0 has to multiplied by the factor 2
     318       751951 :                   e_diff = e_0(ip) - e_old
     319       751951 :                   e_0(ip) = e_0(ip) + e_diff
     320              :                END IF
     321              :             ELSE
     322       393216 :                IF (yval <= 0.0_dp) THEN
     323        10559 :                   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        10559 :                                                sx, R, gamma, grad_deriv)
     327              :                   ! VERY UGLY HACK e_0 has to multiplied by the factor 2
     328        10559 :                   e_diff = e_0(ip) - e_old
     329        10559 :                   e_0(ip) = e_0(ip) + e_diff
     330              :                ELSE
     331       382657 :                   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       382657 :                                               sx, R, gamma, grad_deriv)
     335              :                   ! VERY UGLY HACK e_0 has to multiplied by the factor 2
     336       382657 :                   e_diff = e_0(ip) - e_old
     337       382657 :                   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           84 :    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          410 :    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           82 :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, &
     373           82 :          e_laplace_rhob, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, e_tau_a, e_tau_b, laplace_rhoa, &
     374           82 :          laplace_rhob, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
     375              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     376              : 
     377           82 :       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           82 :                           rho_cutoff=epsilon_rho)
     383           82 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     384              : 
     385           82 :       dummy => rhoa
     386              : 
     387           82 :       e_0 => dummy
     388           82 :       e_rhoa => dummy
     389           82 :       e_rhob => dummy
     390           82 :       e_ndrhoa => dummy
     391           82 :       e_ndrhob => dummy
     392           82 :       e_tau_a => dummy
     393           82 :       e_tau_b => dummy
     394           82 :       e_laplace_rhoa => dummy
     395           82 :       e_laplace_rhob => dummy
     396              : 
     397           82 :       IF (grad_deriv >= 0) THEN
     398              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     399           82 :                                          allocate_deriv=.TRUE.)
     400           82 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     401              :       END IF
     402           82 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     403              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
     404           54 :                                          allocate_deriv=.TRUE.)
     405           54 :          CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
     406              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
     407           54 :                                          allocate_deriv=.TRUE.)
     408           54 :          CALL xc_derivative_get(deriv, deriv_data=e_rhob)
     409              : 
     410              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
     411           54 :                                          allocate_deriv=.TRUE.)
     412           54 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
     413              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
     414           54 :                                          allocate_deriv=.TRUE.)
     415           54 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
     416              : 
     417              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], &
     418           54 :                                          allocate_deriv=.TRUE.)
     419           54 :          CALL xc_derivative_get(deriv, deriv_data=e_tau_a)
     420              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], &
     421           54 :                                          allocate_deriv=.TRUE.)
     422           54 :          CALL xc_derivative_get(deriv, deriv_data=e_tau_b)
     423              : 
     424              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa], &
     425           54 :                                          allocate_deriv=.TRUE.)
     426           54 :          CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa)
     427              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob], &
     428           54 :                                          allocate_deriv=.TRUE.)
     429           54 :          CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob)
     430              :       END IF
     431           82 :       IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
     432            0 :          CPABORT("derivatives bigger than 1 not implemented")
     433              :       END IF
     434              : 
     435           82 :       CALL section_vals_val_get(br_params, "scale_x", r_val=sx)
     436           82 :       CALL section_vals_val_get(br_params, "CUTOFF_RADIUS", r_val=R)
     437           82 :       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           82 : !$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           82 :       CALL timestop(handle)
     462           82 :    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          164 :    SUBROUTINE xbecke_roussel_lsd_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
     490          164 :                                       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      7472250 :          my_rho = MAX(rho(ip), 0.0_dp)
     509      7472250 :          IF (my_rho > epsilon_rho) THEN
     510      7472244 :             my_ndrho = MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp)
     511      7472244 :             my_tau = 1.0_dp*MAX(EPSILON(0.0_dp)*1.e4_dp, tau(ip))
     512      7472244 :             my_laplace_rho = 1.0_dp*laplace_rho(ip)
     513              : 
     514      7472244 :             t1 = pi**(0.1e1_dp/0.3e1_dp)
     515      7472244 :             t2 = t1**2
     516      7472244 :             t3 = my_rho**(0.1e1_dp/0.3e1_dp)
     517      7472244 :             t4 = t3**2
     518      7472244 :             t5 = t4*my_rho
     519      7472244 :             t8 = my_ndrho**2
     520      7472244 :             t9 = 0.1e1_dp/my_rho
     521              :             ! *** CP2K defines tau in a different way as compared to Becke !!!
     522      7472244 :             t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
     523      7472244 :             t16 = 0.1e1_dp/t15
     524      7472244 :             yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
     525      7472244 :             IF (R == 0.0_dp) THEN
     526      2187000 :                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        48750 :                                         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      2138250 :                                        sx, gamma, grad_deriv)
     534              :                END IF
     535              :             ELSE
     536      5285244 :                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       117314 :                                                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      5167930 :                                               sx, R, gamma, grad_deriv)
     544              :                END IF
     545              :             END IF
     546              :          END IF
     547              :       END DO
     548              : 
     549              : !$OMP     END DO
     550              : 
     551          164 :    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       187551 :    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       187551 :       IF (grad_deriv >= 0) THEN
     592       187551 :          t1 = pi**(0.1e1_dp/0.3e1_dp)
     593       187551 :          t2 = t1**2
     594       187551 :          t3 = rho**(0.1e1_dp/0.3e1_dp)
     595       187551 :          t4 = t3**2
     596       187551 :          t5 = t4*rho
     597       187551 :          t9 = ndrho**2
     598       187551 :          t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9/rho/0.4e1_dp)/0.3e1_dp
     599       187551 :          t17 = 0.1e1_dp/t16
     600       187551 :          yval = 0.2e1_dp/0.3e1_dp*t2*t5*t17
     601       187551 :          t19 = t3*rho
     602       187551 :          t20 = 0.3141592654e1_dp**(0.1e1_dp/0.3e1_dp)
     603       187551 :          t21 = t19*t20
     604       187551 :          t22 = br_a1*t2
     605       187551 :          t23 = t5*t17
     606       187551 :          t26 = 0.2e1_dp/0.3e1_dp*t22*t23 + br_a2
     607       187551 :          t27 = ATAN(t26)
     608       187551 :          t28 = -t27 + br_a3
     609       187551 :          t29 = br_c1*t2
     610       187551 :          t32 = t1*pi
     611       187551 :          t33 = br_c2*t32
     612       187551 :          t34 = rho**2
     613       187551 :          t35 = t34*rho
     614       187551 :          t36 = t3*t35
     615       187551 :          t37 = t16**2
     616       187551 :          t38 = 0.1e1_dp/t37
     617       187551 :          t39 = t36*t38
     618       187551 :          t42 = pi**2
     619       187551 :          t43 = br_c3*t42
     620       187551 :          t44 = t34**2
     621       187551 :          t45 = t44*rho
     622       187551 :          t47 = 0.1e1_dp/t37/t16
     623       187551 :          t48 = t45*t47
     624       187551 :          t51 = t2*t42
     625       187551 :          t52 = br_c4*t51
     626       187551 :          t53 = t44*t34
     627       187551 :          t54 = t4*t53
     628       187551 :          t55 = t37**2
     629       187551 :          t56 = 0.1e1_dp/t55
     630       187551 :          t57 = t54*t56
     631       187551 :          t61 = t1*t42*pi
     632       187551 :          t62 = br_c5*t61
     633       187551 :          t63 = t44**2
     634       187551 :          t64 = t3*t63
     635       187551 :          t66 = 0.1e1_dp/t55/t16
     636       187551 :          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       187551 :                /0.243e3_dp*t62*t67
     640       187551 :          t71 = t28*t70
     641       187551 :          t72 = br_b1*t2
     642       187551 :          t75 = br_b2*t32
     643       187551 :          t78 = br_b3*t42
     644       187551 :          t81 = br_b4*t51
     645       187551 :          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       187551 :                /0.243e3_dp*t84*t67
     649       187551 :          t88 = 0.1e1_dp/t87
     650       187551 :          t89 = t71*t88
     651       187551 :          t91 = EXP(t89/0.3e1_dp)
     652       187551 :          t92 = t21*t91
     653       187551 :          t93 = 0.1e1_dp/t28
     654       187551 :          t94 = 0.1e1_dp/t70
     655       187551 :          t95 = t93*t94
     656       187551 :          t96 = EXP(-t89)
     657       187551 :          t97 = t88*t96
     658       187551 :          t100 = 0.1e1_dp - t96 - t71*t97/0.2e1_dp
     659       187551 :          t101 = t87*t100
     660       187551 :          t102 = t95*t101
     661       187551 :          e_0 = e_0 + (-t92*t102)*sx
     662              :       END IF
     663       187551 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     664       187551 :          t108 = t4*t17
     665       187551 :          t111 = 0.1e1_dp/t3
     666       187551 :          t113 = t38*gamma
     667       187551 :          t114 = t113*t9
     668       187551 :          t117 = 0.10e2_dp/0.9e1_dp*t22*t108 + t22*t111*t114/0.18e2_dp
     669       187551 :          t118 = t26**2
     670       187551 :          t120 = 0.1e1_dp/(0.1e1_dp + t118)
     671       187551 :          t121 = t117*t120
     672       187551 :          t122 = t70*t88
     673       187551 :          t129 = t3*t34
     674       187551 :          t130 = t129*t38
     675       187551 :          t134 = t47*gamma
     676       187551 :          t135 = t134*t9
     677       187551 :          t138 = t44*t47
     678       187551 :          t142 = t56*gamma
     679       187551 :          t143 = t142*t9
     680       187551 :          t146 = t4*t45
     681       187551 :          t147 = t146*t56
     682       187551 :          t150 = t4*t44
     683       187551 :          t152 = t66*gamma
     684       187551 :          t153 = t152*t9
     685       187551 :          t157 = t3*t44*t35
     686       187551 :          t158 = t157*t66
     687       187551 :          t161 = t3*t53
     688       187551 :          t164 = 0.1e1_dp/t55/t37
     689       187551 :          t165 = t164*gamma
     690       187551 :          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       187551 :                 *t166
     697       187551 :          t170 = t28*t169
     698       187551 :          t172 = t87**2
     699       187551 :          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       187551 :                 *t166
     706       187551 :          t202 = -t121*t122 + t170*t88 - t71*t173*t199
     707       187551 :          t207 = t28**2
     708       187551 :          t208 = 0.1e1_dp/t207
     709       187551 :          t209 = t91*t208
     710       187551 :          t217 = t21*t91*t93
     711       187551 :          t218 = t70**2
     712       187551 :          t220 = 0.1e1_dp/t218*t87
     713       187551 :          t227 = -t202
     714       187551 :          t229 = t122*t96
     715       187551 :          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       187551 :                            *t199/0.2e1_dp - t71*t88*t227*t96/0.2e1_dp))*sx
     721       187551 :          t246 = t4*t38
     722       187551 :          t249 = t120*t70
     723       187551 :          t252 = t22*t246*gamma*ndrho*t249*t88
     724       187551 :          t255 = t113*ndrho
     725       187551 :          t259 = t134*ndrho
     726       187551 :          t263 = t142*ndrho
     727       187551 :          t267 = t152*ndrho
     728       187551 :          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       187551 :                 *t267 - 0.80e2_dp/0.729e3_dp*t62*t157*t271
     732       187551 :          t275 = t28*t274
     733       187551 :          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       187551 :                 *t267 - 0.80e2_dp/0.729e3_dp*t84*t157*t271
     737       187551 :          t295 = t71*t173*t293
     738       187551 :          t304 = t208*t94*t87
     739       187551 :          t307 = t100*br_a1*t2
     740       187551 :          t308 = ndrho*t120
     741       187551 :          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       187551 :                                 /0.2e1_dp))*sx
     748       187551 :          t340 = t5*t38
     749       187551 :          t341 = t22*t340
     750       187551 :          t342 = gamma*t120
     751       187551 :          t344 = t341*t342*t122
     752       187551 :          t346 = t340*gamma
     753       187551 :          t349 = t36*t47
     754       187551 :          t350 = t349*gamma
     755       187551 :          t353 = t45*t56
     756       187551 :          t354 = t353*gamma
     757       187551 :          t357 = t54*t66
     758       187551 :          t358 = t357*gamma
     759       187551 :          t361 = t64*t164
     760       187551 :          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       187551 :                 /0.729e3_dp*t62*t362
     764       187551 :          t366 = t28*t365
     765       187551 :          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       187551 :                 /0.729e3_dp*t84*t362
     769       187551 :          t381 = t71*t173*t379
     770       187551 :          t387 = t35*t20
     771       187551 :          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       187551 :                                     *t96/0.2e1_dp))*sx
     778       187551 :          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       187551 :                 *t62*t361
     782       187551 :          t435 = t28*t434
     783       187551 :          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       187551 :                 *t84*t361
     787       187551 :          t450 = t71*t173*t448
     788       187551 :          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       187551 :                                                 + t71*t234*t448/0.2e1_dp - t71*t88*t471*t96/0.2e1_dp))*sx
     794              :       END IF
     795       187551 :    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      2890201 :    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      2890201 :       IF (grad_deriv >= 0) THEN
     836      2890201 :          t1 = pi**(0.1e1_dp/0.3e1_dp)
     837      2890201 :          t2 = t1**2
     838      2890201 :          t3 = rho**(0.1e1_dp/0.3e1_dp)
     839      2890201 :          t4 = t3**2
     840      2890201 :          t5 = t4*rho
     841      2890201 :          t6 = t2*t5
     842      2890201 :          t9 = ndrho**2
     843      2890201 :          t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9/rho/0.4e1_dp)/0.3e1_dp
     844      2890201 :          t17 = 0.1e1_dp/t16
     845      2890201 :          yval = 0.2e1_dp/0.3e1_dp*t6*t17
     846      2890201 :          t19 = t3*rho
     847      2890201 :          t20 = 0.3141592654e1_dp**(0.1e1_dp/0.3e1_dp)
     848      2890201 :          t21 = t19*t20
     849      2890201 :          t22 = 0.1e1_dp/br_BB
     850      2890201 :          t23 = 0.1e1_dp/t2
     851      2890201 :          t24 = t22*t23
     852      2890201 :          t25 = 0.1e1_dp/t5
     853      2890201 :          t26 = t25*t16
     854      2890201 :          t29 = br_BB**2
     855      2890201 :          t30 = t1*pi
     856      2890201 :          t31 = t29*t30
     857      2890201 :          t32 = rho**2
     858      2890201 :          t33 = t32*rho
     859      2890201 :          t34 = t3*t33
     860      2890201 :          t35 = t16**2
     861      2890201 :          t36 = 0.1e1_dp/t35
     862      2890201 :          t37 = t34*t36
     863      2890201 :          t41 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t31*t37)
     864      2890201 :          t42 = t41*t22
     865      2890201 :          t43 = t23*t25
     866              :          t47 = 0.1500000000000000e1_dp*t24*t26 + 0.3e1_dp/0.2e1_dp*t42*t43 &
     867      2890201 :                *t16
     868      2890201 :          t48 = LOG(t47)
     869      2890201 :          t49 = t48 + 0.2e1_dp
     870      2890201 :          t50 = br_d1*t2
     871      2890201 :          t51 = t5*t17
     872      2890201 :          t54 = br_d2*t30
     873      2890201 :          t57 = pi**2
     874      2890201 :          t58 = br_d3*t57
     875      2890201 :          t59 = t32**2
     876      2890201 :          t60 = t59*rho
     877      2890201 :          t62 = 0.1e1_dp/t35/t16
     878      2890201 :          t63 = t60*t62
     879      2890201 :          t66 = t2*t57
     880      2890201 :          t67 = br_d4*t66
     881      2890201 :          t68 = t59*t32
     882      2890201 :          t69 = t4*t68
     883      2890201 :          t70 = t35**2
     884      2890201 :          t71 = 0.1e1_dp/t70
     885      2890201 :          t72 = t69*t71
     886      2890201 :          t76 = t1*t57*pi
     887      2890201 :          t77 = br_d5*t76
     888      2890201 :          t78 = t59**2
     889      2890201 :          t79 = t3*t78
     890      2890201 :          t81 = 0.1e1_dp/t70/t16
     891      2890201 :          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      2890201 :                /0.243e3_dp*t77*t82
     895      2890201 :          t86 = t49*t85
     896      2890201 :          t87 = br_e1*t2
     897      2890201 :          t90 = br_e2*t30
     898      2890201 :          t93 = br_e3*t57
     899      2890201 :          t96 = br_e4*t66
     900      2890201 :          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      2890201 :                 /0.243e3_dp*t99*t82
     904      2890201 :          t103 = 0.1e1_dp/t102
     905      2890201 :          t104 = t86*t103
     906      2890201 :          t106 = EXP(t104/0.3e1_dp)
     907      2890201 :          t107 = t21*t106
     908      2890201 :          t108 = 0.1e1_dp/t49
     909      2890201 :          t109 = 0.1e1_dp/t85
     910      2890201 :          t110 = t108*t109
     911      2890201 :          t111 = EXP(-t104)
     912      2890201 :          t112 = t103*t111
     913      2890201 :          t115 = 0.1e1_dp - t111 - t86*t112/0.2e1_dp
     914      2890201 :          t117 = t110*t102*t115
     915      2890201 :          e_0 = e_0 + (-t107*t117)*sx
     916              :       END IF
     917      2890201 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     918      2890201 :          t124 = 0.1e1_dp/t4/t32
     919      2890201 :          t131 = 0.1e1_dp/t4/t33*gamma*t9
     920      2890201 :          t134 = 0.1e1_dp/t41
     921      2890201 :          t137 = t3*t32
     922      2890201 :          t138 = t137*t36
     923      2890201 :          t142 = t62*gamma
     924      2890201 :          t143 = t142*t9
     925      2890201 :          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      2890201 :                 t42*t23*t124*t16 - t154*t131/0.8e1_dp
     930      2890201 :          t158 = 0.1e1_dp/t47
     931      2890201 :          t159 = t157*t158
     932      2890201 :          t160 = t85*t103
     933      2890201 :          t162 = t4*t17
     934      2890201 :          t165 = 0.1e1_dp/t3
     935      2890201 :          t167 = t36*gamma
     936      2890201 :          t168 = t167*t9
     937      2890201 :          t176 = t59*t62
     938      2890201 :          t180 = t71*gamma
     939      2890201 :          t181 = t180*t9
     940      2890201 :          t184 = t4*t60
     941      2890201 :          t185 = t184*t71
     942      2890201 :          t188 = t4*t59
     943      2890201 :          t190 = t81*gamma
     944      2890201 :          t191 = t190*t9
     945      2890201 :          t195 = t3*t59*t33
     946      2890201 :          t196 = t195*t81
     947      2890201 :          t199 = t3*t68
     948      2890201 :          t202 = 0.1e1_dp/t70/t35
     949      2890201 :          t203 = t202*gamma
     950      2890201 :          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      2890201 :                 *t204
     957      2890201 :          t208 = t49*t207
     958      2890201 :          t210 = t102**2
     959      2890201 :          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      2890201 :                 *t204
     966      2890201 :          t240 = t159*t160 + t208*t103 - t86*t211*t237
     967      2890201 :          t245 = t49**2
     968      2890201 :          t248 = t21*t106/t245
     969      2890201 :          t249 = t109*t102
     970      2890201 :          t255 = t21*t106*t108
     971      2890201 :          t256 = t85**2
     972      2890201 :          t258 = 0.1e1_dp/t256*t102
     973      2890201 :          t265 = -t240
     974      2890201 :          t267 = t160*t111
     975      2890201 :          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      2890201 :                                                                             *t237/0.2e1_dp - t86*t103*t265*t111/0.2e1_dp))*sx
     981      2890201 :          t285 = t124*gamma*ndrho
     982      2890201 :          t288 = t134*br_BB
     983      2890201 :          t289 = t288*t2
     984              :          t297 = 0.2500000000000000e0_dp*t24*t285 - t289*t4*t36*gamma* &
     985      2890201 :                 ndrho/0.9e1_dp + t154*t285/0.4e1_dp
     986      2890201 :          t298 = t297*t158
     987      2890201 :          t301 = t167*ndrho
     988      2890201 :          t305 = t142*ndrho
     989      2890201 :          t309 = t180*ndrho
     990      2890201 :          t313 = t190*ndrho
     991      2890201 :          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      2890201 :                 *t313 - 0.80e2_dp/0.729e3_dp*t77*t195*t317
     995      2890201 :          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      2890201 :                 *t313 - 0.80e2_dp/0.729e3_dp*t99*t195*t317
     999      2890201 :          t341 = t298*t160 + t321*t103 - t86*t211*t338
    1000      2890201 :          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      2890201 :                                                 t111/0.2e1_dp))*sx
    1006      2890201 :          t376 = t5*t36
    1007      2890201 :          t377 = t376*gamma
    1008              :          t382 = -0.1000000000e1_dp*t24*t25*gamma + 0.4e1_dp/0.9e1_dp*t289*t377 &
    1009      2890201 :                 - t42*t43*gamma
    1010      2890201 :          t383 = t382*t158
    1011      2890201 :          t387 = t34*t62
    1012      2890201 :          t388 = t387*gamma
    1013      2890201 :          t391 = t60*t71
    1014      2890201 :          t392 = t391*gamma
    1015      2890201 :          t395 = t69*t81
    1016      2890201 :          t396 = t395*gamma
    1017      2890201 :          t399 = t79*t202
    1018      2890201 :          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      2890201 :                 /0.729e3_dp*t77*t400
    1022      2890201 :          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      2890201 :                 /0.729e3_dp*t99*t400
    1026      2890201 :          t419 = t383*t160 + t404*t103 - t86*t211*t416
    1027      2890201 :          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      2890201 :                                           /0.2e1_dp))*sx
    1033              :          t458 = 0.2500000000000000e0_dp*t24*t25 - t288*t6*t36/0.9e1_dp + &
    1034      2890201 :                 t42*t43/0.4e1_dp
    1035      2890201 :          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      2890201 :                 *t77*t399
    1039      2890201 :          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      2890201 :                 *t99*t399
    1043      2890201 :          t487 = t459*t160 + t472*t103 - t86*t211*t484
    1044      2890201 :          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      2890201 :                                                           /0.2e1_dp))*sx
    1050              :       END IF
    1051              : 
    1052      2890201 :    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       386143 :    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       386143 :       t1 = 8**(0.1e1_dp/0.3e1_dp)
    1088       386143 :       t2 = t1**2
    1089       386143 :       t3 = pi**(0.1e1_dp/0.3e1_dp)
    1090       386143 :       t4 = t3**2
    1091       386143 :       t6 = rho**(0.1e1_dp/0.3e1_dp)
    1092       386143 :       t7 = t6**2
    1093       386143 :       t11 = ndrho**2
    1094       386143 :       t12 = 0.1e1_dp/rho
    1095       386143 :       t18 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t11*t12/0.4e1_dp)/0.3e1_dp
    1096       386143 :       t20 = t7*rho/t18
    1097       386143 :       t24 = ATAN(0.2e1_dp/0.3e1_dp*br_a1*t4*t20 + br_a2)
    1098       386143 :       t25 = -t24 + br_a3
    1099       386143 :       t26 = t25**2
    1100       386143 :       t31 = t3*pi
    1101       386143 :       t33 = rho**2
    1102       386143 :       t36 = t18**2
    1103       386143 :       t38 = t6*t33*rho/t36
    1104       386143 :       t41 = pi**2
    1105       386143 :       t43 = t33**2
    1106       386143 :       t47 = t43*rho/t36/t18
    1107       386143 :       t50 = t4*t41
    1108       386143 :       t54 = t36**2
    1109       386143 :       t56 = t7*t43*t33/t54
    1110       386143 :       t60 = t3*t41*pi
    1111       386143 :       t62 = t43**2
    1112       386143 :       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       386143 :             *br_c4*t50*t56 + 0.32e2_dp/0.243e3_dp*br_c5*t60*t66
    1116       386143 :       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       386143 :             *br_b4*t50*t56 + 0.32e2_dp/0.243e3_dp*br_b5*t60*t66
    1120       386143 :       t89 = t88**2
    1121       386143 :       t96 = EXP(-t25*t69/t88)
    1122              :       t101 = (t26*t25*t70*t69/t89/t88*t96/0.3141592654e1_dp* &
    1123       386143 :               t12)**(0.1e1_dp/0.3e1_dp)
    1124       386143 :       brval = REAL(t2, KIND=dp)*t101/0.8e1_dp
    1125              : 
    1126       386143 :       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       368598 :                                              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        17545 :                                               sx, R, gamma, grad_deriv)
    1134              :       END IF
    1135              : 
    1136       386143 :    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       368598 :    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       368598 :       IF (grad_deriv >= 0) THEN
    1183       368598 :          t1 = pi**(0.1e1_dp/0.3e1_dp)
    1184       368598 :          t2 = t1**2
    1185       368598 :          t3 = br_a1*t2
    1186       368598 :          t4 = rho**(0.1e1_dp/0.3e1_dp)
    1187       368598 :          t5 = t4**2
    1188       368598 :          t6 = t5*rho
    1189       368598 :          t9 = ndrho**2
    1190       368598 :          t10 = 0.1e1_dp/rho
    1191       368598 :          t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9*t10/0.4e1_dp)/0.3e1_dp
    1192       368598 :          t17 = 0.1e1_dp/t16
    1193       368598 :          t18 = t6*t17
    1194       368598 :          t21 = 0.2e1_dp/0.3e1_dp*t3*t18 + br_a2
    1195       368598 :          t22 = ATAN(t21)
    1196       368598 :          t23 = -t22 + br_a3
    1197       368598 :          t24 = br_c1*t2
    1198       368598 :          t27 = t1*pi
    1199       368598 :          t28 = br_c2*t27
    1200       368598 :          t29 = rho**2
    1201       368598 :          t30 = t29*rho
    1202       368598 :          t31 = t4*t30
    1203       368598 :          t32 = t16**2
    1204       368598 :          t33 = 0.1e1_dp/t32
    1205       368598 :          t34 = t31*t33
    1206       368598 :          t37 = pi**2
    1207       368598 :          t38 = br_c3*t37
    1208       368598 :          t39 = t29**2
    1209       368598 :          t40 = t39*rho
    1210       368598 :          t42 = 0.1e1_dp/t32/t16
    1211       368598 :          t43 = t40*t42
    1212       368598 :          t46 = t2*t37
    1213       368598 :          t47 = br_c4*t46
    1214       368598 :          t48 = t39*t29
    1215       368598 :          t49 = t5*t48
    1216       368598 :          t50 = t32**2
    1217       368598 :          t51 = 0.1e1_dp/t50
    1218       368598 :          t52 = t49*t51
    1219       368598 :          t56 = t1*t37*pi
    1220       368598 :          t57 = br_c5*t56
    1221       368598 :          t58 = t39**2
    1222       368598 :          t59 = t4*t58
    1223       368598 :          t61 = 0.1e1_dp/t50/t16
    1224       368598 :          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       368598 :                /0.243e3_dp*t57*t62
    1228       368598 :          t66 = t23*t65
    1229       368598 :          t67 = br_b1*t2
    1230       368598 :          t70 = br_b2*t27
    1231       368598 :          t73 = br_b3*t37
    1232       368598 :          t76 = br_b4*t46
    1233       368598 :          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       368598 :                /0.243e3_dp*t79*t62
    1237       368598 :          t83 = 0.1e1_dp/t82
    1238       368598 :          t84 = t66*t83
    1239       368598 :          t85 = 8**(0.1e1_dp/0.3e1_dp)
    1240       368598 :          t86 = t23**2
    1241       368598 :          t87 = t86*t23
    1242       368598 :          t88 = t65**2
    1243       368598 :          t89 = t88*t65
    1244       368598 :          t90 = t87*t89
    1245       368598 :          t91 = t82**2
    1246       368598 :          t93 = 0.1e1_dp/t91/t82
    1247       368598 :          t94 = t90*t93
    1248       368598 :          t95 = EXP(-t84)
    1249       368598 :          t96 = 0.1e1_dp/0.3141592654e1_dp
    1250       368598 :          t97 = t95*t96
    1251       368598 :          t99 = t94*t97*t10
    1252       368598 :          t100 = t99**(0.1e1_dp/0.3e1_dp)
    1253       368598 :          t101 = 0.1e1_dp/t100
    1254       368598 :          t102 = REAL(t85, KIND=dp)*t101
    1255       368598 :          t103 = t102*R
    1256       368598 :          t104 = t84*t103
    1257       368598 :          t106 = EXP(t84 - t104)
    1258       368598 :          t108 = t106*t23
    1259       368598 :          t109 = t65*t83
    1260       368598 :          t110 = t108*t109
    1261       368598 :          t114 = t83*REAL(t85, KIND=dp)*t101*R
    1262       368598 :          t117 = EXP(-t84 - t104)
    1263       368598 :          t119 = t117*t23
    1264       368598 :          t120 = t119*t109
    1265              :          t123 = -0.2e1_dp*t106 + t110 - t108*t65*t114 + 0.2e1_dp*t117 + t120 &
    1266       368598 :                 + t119*t65*t114
    1267       368598 :          t124 = rho*t123
    1268       368598 :          e_0 = e_0 + (t124*t102/0.8e1_dp)*sx
    1269              :       END IF
    1270       368598 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1271       194202 :          t129 = t5*t17
    1272       194202 :          t132 = 0.1e1_dp/t4
    1273       194202 :          t134 = t33*gamma
    1274       194202 :          t135 = t134*t9
    1275       194202 :          t138 = 0.10e2_dp/0.9e1_dp*t3*t129 + t3*t132*t135/0.18e2_dp
    1276       194202 :          t139 = t21**2
    1277       194202 :          t141 = 0.1e1_dp/(0.1e1_dp + t139)
    1278       194202 :          t142 = t138*t141
    1279       194202 :          t143 = t142*t109
    1280       194202 :          t149 = t4*t29
    1281       194202 :          t150 = t149*t33
    1282       194202 :          t153 = t4*rho
    1283       194202 :          t155 = t42*gamma
    1284       194202 :          t156 = t155*t9
    1285       194202 :          t159 = t39*t42
    1286       194202 :          t163 = t51*gamma
    1287       194202 :          t164 = t163*t9
    1288       194202 :          t167 = t5*t40
    1289       194202 :          t168 = t167*t51
    1290       194202 :          t171 = t5*t39
    1291       194202 :          t173 = t61*gamma
    1292       194202 :          t174 = t173*t9
    1293       194202 :          t178 = t4*t39*t30
    1294       194202 :          t179 = t178*t61
    1295       194202 :          t182 = t4*t48
    1296       194202 :          t185 = 0.1e1_dp/t50/t32
    1297       194202 :          t186 = t185*gamma
    1298       194202 :          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       194202 :                 t182*t187
    1305       194202 :          t192 = t23*t190*t83
    1306       194202 :          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       194202 :                 t182*t187
    1313       194202 :          t221 = t66*t193*t219
    1314       194202 :          t223 = t142*t65*t114
    1315       194202 :          t224 = t192*t103
    1316       194202 :          t225 = t66*t193
    1317       194202 :          t227 = t102*R*t219
    1318       194202 :          t228 = t225*t227
    1319       194202 :          t231 = REAL(t85, KIND=dp)/t100/t99
    1320       194202 :          t232 = t86*t89
    1321       194202 :          t233 = t93*t95
    1322       194202 :          t235 = t96*t10
    1323       194202 :          t240 = t87*t88*t93
    1324       194202 :          t245 = t91**2
    1325       194202 :          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       194202 :                                                           t221)*t95*t235 - t94*t97/t29
    1329       194202 :          t261 = t231*R*t259
    1330       194202 :          t263 = t84*t261/0.3e1_dp
    1331       194202 :          t265 = (-t143 + t192 - t221 + t223 - t224 + t228 + t263)*t106
    1332       194202 :          t268 = t106*t138
    1333       194202 :          t269 = t141*t65
    1334       194202 :          t270 = t269*t83
    1335       194202 :          t272 = t190*t83
    1336       194202 :          t274 = t65*t193
    1337       194202 :          t275 = t274*t219
    1338       194202 :          t283 = t108*t274
    1339       194202 :          t288 = (t143 - t192 + t221 + t223 - t224 + t228 + t263)*t117
    1340       194202 :          t291 = t117*t138
    1341       194202 :          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       194202 :                 /0.3e1_dp
    1348              :          e_rho = e_rho + (t123*REAL(t85, KIND=dp)*t101/0.8e1_dp + rho*t305*t102/0.8e1_dp &
    1349       194202 :                           - t124*t231*t259/0.24e2_dp)*sx
    1350       194202 :          t312 = t5*t33
    1351       194202 :          t314 = gamma*ndrho
    1352       194202 :          t315 = t314*t270
    1353       194202 :          t317 = t3*t312*t315/0.9e1_dp
    1354       194202 :          t319 = t134*ndrho
    1355       194202 :          t323 = t155*ndrho
    1356       194202 :          t327 = t163*ndrho
    1357       194202 :          t331 = t173*ndrho
    1358       194202 :          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       194202 :                 *t331 - 0.80e2_dp/0.729e3_dp*t57*t178*t335
    1362       194202 :          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       194202 :                 *t331 - 0.80e2_dp/0.729e3_dp*t79*t178*t335
    1366       194202 :          t358 = t66*t193*t356
    1367       194202 :          t359 = t3*t5
    1368       194202 :          t361 = t270*t103
    1369       194202 :          t363 = t359*t319*t361/0.9e1_dp
    1370       194202 :          t364 = t340*t103
    1371       194202 :          t366 = t102*R*t356
    1372       194202 :          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       194202 :                 t94*(-t317 - t340 + t358)*t95*t235
    1376       194202 :          t390 = t231*R*t388
    1377       194202 :          t392 = t84*t390/0.3e1_dp
    1378       194202 :          t394 = (t317 + t340 - t358 - t363 - t364 + t367 + t392)*t106
    1379       194202 :          t397 = t106*br_a1
    1380       194202 :          t399 = t2*t5*t33
    1381       194202 :          t403 = t338*t83
    1382       194202 :          t405 = t274*t356
    1383       194202 :          t409 = t397*t2
    1384       194202 :          t410 = t312*gamma
    1385       194202 :          t414 = ndrho*t141*t65*t114
    1386       194202 :          t423 = (-t317 - t340 + t358 - t363 - t364 + t367 + t392)*t117
    1387       194202 :          t426 = t117*br_a1
    1388       194202 :          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       194202 :                 + t119*t338*t114 - t301*t366 - t120*t390/0.3e1_dp
    1395       194202 :          e_ndrho = e_ndrho + (rho*t443*t102/0.8e1_dp - t124*t231*t388/0.24e2_dp)*sx
    1396       194202 :          t450 = t6*t33
    1397       194202 :          t455 = 0.4e1_dp/0.9e1_dp*t3*t450*gamma*t141*t109
    1398       194202 :          t456 = t450*gamma
    1399       194202 :          t459 = t31*t42
    1400       194202 :          t460 = t459*gamma
    1401       194202 :          t463 = t40*t51
    1402       194202 :          t464 = t463*gamma
    1403       194202 :          t467 = t49*t61
    1404       194202 :          t468 = t467*gamma
    1405       194202 :          t471 = t59*t185
    1406       194202 :          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       194202 :                 /0.729e3_dp*t57*t472
    1410       194202 :          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       194202 :                 /0.729e3_dp*t79*t472
    1414       194202 :          t490 = t66*t193*t488
    1415       194202 :          t493 = 0.4e1_dp/0.9e1_dp*t3*t456*t361
    1416       194202 :          t494 = t477*t103
    1417       194202 :          t496 = t102*R*t488
    1418       194202 :          t497 = t225*t496
    1419       194202 :          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       194202 :                                                                  t477 + t490)*t95*t235
    1423       194202 :          t518 = t231*R*t516
    1424       194202 :          t520 = t84*t518/0.3e1_dp
    1425       194202 :          t522 = (-t455 + t477 - t490 + t493 - t494 + t497 + t520)*t106
    1426       194202 :          t525 = t2*t6
    1427       194202 :          t526 = t397*t525
    1428       194202 :          t527 = t134*t270
    1429       194202 :          t530 = t475*t83
    1430       194202 :          t532 = t274*t488
    1431       194202 :          t545 = (t455 - t477 + t490 + t493 - t494 + t497 + t520)*t117
    1432       194202 :          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       194202 :                 /0.3e1_dp
    1440       194202 :          e_tau = e_tau + (rho*t563*t102/0.8e1_dp - t124*t231*t516/0.24e2_dp)*sx
    1441       194202 :          t572 = t33*t141*t109
    1442       194202 :          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       194202 :                 *t57*t471
    1446       194202 :          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       194202 :                 *t79*t471
    1450       194202 :          t600 = t66*t193*t598
    1451       194202 :          t605 = t3*t450*t141*t109*t103/0.9e1_dp
    1452       194202 :          t606 = t587*t103
    1453       194202 :          t608 = t102*R*t598
    1454       194202 :          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       194202 :                                                                  - t587 + t600)*t95*t235
    1458       194202 :          t630 = t231*R*t628
    1459       194202 :          t632 = t84*t630/0.3e1_dp
    1460       194202 :          t634 = (t574 + t587 - t600 - t605 - t606 + t609 + t632)*t106
    1461       194202 :          t639 = t585*t83
    1462       194202 :          t641 = t274*t598
    1463       194202 :          t645 = t525*t33
    1464       194202 :          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       194202 :                 *t114 - t301*t608 - t120*t630/0.3e1_dp
    1471       194202 :          e_laplace_rho = e_laplace_rho + (rho*t672*t102/0.8e1_dp - t124*t231*t628/0.24e2_dp)*sx
    1472              :       END IF
    1473              : 
    1474       368598 :    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        17545 :    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        17545 :       IF (grad_deriv >= 0) THEN
    1521        17545 :          t1 = pi**(0.1e1_dp/0.3e1_dp)
    1522        17545 :          t2 = t1**2
    1523        17545 :          t3 = br_a1*t2
    1524        17545 :          t4 = rho**(0.1e1_dp/0.3e1_dp)
    1525        17545 :          t5 = t4**2
    1526        17545 :          t6 = t5*rho
    1527        17545 :          t9 = ndrho**2
    1528        17545 :          t10 = 0.1e1_dp/rho
    1529        17545 :          t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9*t10/0.4e1_dp)/0.3e1_dp
    1530        17545 :          t17 = 0.1e1_dp/t16
    1531        17545 :          t18 = t6*t17
    1532        17545 :          t21 = 0.2e1_dp/0.3e1_dp*t3*t18 + br_a2
    1533        17545 :          t22 = ATAN(t21)
    1534        17545 :          t23 = -t22 + br_a3
    1535        17545 :          t24 = br_c1*t2
    1536        17545 :          t27 = t1*pi
    1537        17545 :          t28 = br_c2*t27
    1538        17545 :          t29 = rho**2
    1539        17545 :          t30 = t29*rho
    1540        17545 :          t31 = t4*t30
    1541        17545 :          t32 = t16**2
    1542        17545 :          t33 = 0.1e1_dp/t32
    1543        17545 :          t34 = t31*t33
    1544        17545 :          t37 = pi**2
    1545        17545 :          t38 = br_c3*t37
    1546        17545 :          t39 = t29**2
    1547        17545 :          t40 = t39*rho
    1548        17545 :          t42 = 0.1e1_dp/t32/t16
    1549        17545 :          t43 = t40*t42
    1550        17545 :          t46 = t2*t37
    1551        17545 :          t47 = br_c4*t46
    1552        17545 :          t48 = t39*t29
    1553        17545 :          t49 = t5*t48
    1554        17545 :          t50 = t32**2
    1555        17545 :          t51 = 0.1e1_dp/t50
    1556        17545 :          t52 = t49*t51
    1557        17545 :          t56 = t1*t37*pi
    1558        17545 :          t57 = br_c5*t56
    1559        17545 :          t58 = t39**2
    1560        17545 :          t59 = t4*t58
    1561        17545 :          t61 = 0.1e1_dp/t50/t16
    1562        17545 :          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        17545 :                /0.243e3_dp*t57*t62
    1566        17545 :          t66 = t23*t65
    1567        17545 :          t67 = br_b1*t2
    1568        17545 :          t70 = br_b2*t27
    1569        17545 :          t73 = br_b3*t37
    1570        17545 :          t76 = br_b4*t46
    1571        17545 :          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        17545 :                /0.243e3_dp*t79*t62
    1575        17545 :          t83 = 0.1e1_dp/t82
    1576        17545 :          t84 = t66*t83
    1577        17545 :          t85 = 8**(0.1e1_dp/0.3e1_dp)
    1578        17545 :          t86 = t23**2
    1579        17545 :          t87 = t86*t23
    1580        17545 :          t88 = t65**2
    1581        17545 :          t89 = t88*t65
    1582        17545 :          t90 = t87*t89
    1583        17545 :          t91 = t82**2
    1584        17545 :          t93 = 0.1e1_dp/t91/t82
    1585        17545 :          t94 = t90*t93
    1586        17545 :          t95 = EXP(-t84)
    1587        17545 :          t96 = 0.1e1_dp/0.3141592654e1_dp
    1588        17545 :          t97 = t95*t96
    1589        17545 :          t99 = t94*t97*t10
    1590        17545 :          t100 = t99**(0.1e1_dp/0.3e1_dp)
    1591        17545 :          t101 = 0.1e1_dp/t100
    1592        17545 :          t102 = REAL(t85, KIND=dp)*t101
    1593        17545 :          t103 = t102*R
    1594        17545 :          t104 = t84*t103
    1595        17545 :          t106 = EXP(0.2e1_dp*t104)
    1596        17545 :          t108 = t106*R
    1597        17545 :          t109 = t108*t23
    1598        17545 :          t110 = t65*t83
    1599        17545 :          t111 = t110*t102
    1600        17545 :          t113 = t106*t23
    1601        17545 :          t115 = t84 + t104
    1602        17545 :          t116 = EXP(t115)
    1603              :          t118 = 0.2e1_dp*t106 - t109*t111 + t113*t110 + 0.2e1_dp + t84 + t104 &
    1604        17545 :                 - 0.4e1_dp*t116
    1605        17545 :          t119 = rho*t118
    1606        17545 :          t121 = EXP(-t115)
    1607        17545 :          t123 = t121*REAL(t85, KIND=dp)*t101
    1608        17545 :          e_0 = e_0 + (t119*t123/0.8e1_dp)*sx
    1609              :       END IF
    1610        17545 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1611         9431 :          t128 = t5*t17
    1612         9431 :          t131 = 0.1e1_dp/t4
    1613         9431 :          t133 = t33*gamma
    1614         9431 :          t134 = t133*t9
    1615         9431 :          t137 = 0.10e2_dp/0.9e1_dp*t3*t128 + t3*t131*t134/0.18e2_dp
    1616         9431 :          t138 = t21**2
    1617         9431 :          t140 = 0.1e1_dp/(0.1e1_dp + t138)
    1618         9431 :          t141 = t137*t140
    1619         9431 :          t143 = t83*REAL(t85, KIND=dp)
    1620         9431 :          t146 = t141*t65*t143*t101*R
    1621         9431 :          t153 = t4*t29
    1622         9431 :          t154 = t153*t33
    1623         9431 :          t157 = t4*rho
    1624         9431 :          t159 = t42*gamma
    1625         9431 :          t160 = t159*t9
    1626         9431 :          t163 = t39*t42
    1627         9431 :          t167 = t51*gamma
    1628         9431 :          t168 = t167*t9
    1629         9431 :          t171 = t5*t40
    1630         9431 :          t172 = t171*t51
    1631         9431 :          t175 = t5*t39
    1632         9431 :          t177 = t61*gamma
    1633         9431 :          t178 = t177*t9
    1634         9431 :          t182 = t4*t39*t30
    1635         9431 :          t183 = t182*t61
    1636         9431 :          t186 = t4*t48
    1637         9431 :          t189 = 0.1e1_dp/t50/t32
    1638         9431 :          t190 = t189*gamma
    1639         9431 :          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         9431 :                 t186*t191
    1646         9431 :          t196 = t23*t194*t83
    1647         9431 :          t197 = t196*t103
    1648         9431 :          t199 = 0.1e1_dp/t91
    1649         9431 :          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         9431 :                 t186*t191
    1656         9431 :          t229 = t200*t102*R*t226
    1657         9431 :          t232 = 0.1e1_dp/t100/t99
    1658         9431 :          t233 = REAL(t85, KIND=dp)*t232
    1659         9431 :          t234 = t86*t89
    1660         9431 :          t235 = t93*t95
    1661         9431 :          t237 = t96*t10
    1662         9431 :          t242 = t87*t88*t93
    1663         9431 :          t247 = t91**2
    1664         9431 :          t249 = t90/t247
    1665         9431 :          t254 = t141*t110
    1666         9431 :          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         9431 :                                                           t256)*t95*t237 - t94*t97/t29
    1670         9431 :          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         9431 :                  *t267)*t106
    1673         9431 :          t272 = R*t23
    1674         9431 :          t277 = t194*t83
    1675         9431 :          t280 = t108*t66
    1676         9431 :          t281 = t199*REAL(t85, KIND=dp)
    1677         9431 :          t292 = t140*t65*t83
    1678         9431 :          t295 = t65*t199
    1679         9431 :          t298 = t267/0.3e1_dp
    1680         9431 :          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         9431 :                 - 0.4e1_dp*t299*t116
    1686         9431 :          t310 = t119*t121
    1687              :          e_rho = e_rho + (t118*t121*t102/0.8e1_dp + rho*t302*t123/0.8e1_dp - t119 &
    1688         9431 :                           *t299*t123/0.8e1_dp - t310*t233*t264/0.24e2_dp)*sx
    1689         9431 :          t314 = t3*t5
    1690         9431 :          t315 = t133*ndrho
    1691         9431 :          t317 = t292*t103
    1692         9431 :          t318 = t314*t315*t317
    1693         9431 :          t324 = t159*ndrho
    1694         9431 :          t328 = t167*ndrho
    1695         9431 :          t332 = t177*ndrho
    1696         9431 :          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         9431 :                 *t332 - 0.80e2_dp/0.729e3_dp*t57*t182*t336
    1700         9431 :          t341 = t23*t339*t83
    1701         9431 :          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         9431 :                 *t332 - 0.80e2_dp/0.729e3_dp*t79*t182*t336
    1705         9431 :          t362 = t200*t102*R*t359
    1706         9431 :          t368 = gamma*ndrho
    1707         9431 :          t369 = t368*t140
    1708         9431 :          t383 = t368*t292
    1709         9431 :          t385 = t3*t5*t33*t383/0.9e1_dp
    1710         9431 :          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         9431 :                 (-t385 - t341 + t387)*t95*t237
    1714         9431 :          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         9431 :                  /0.3e1_dp*t395)*t106
    1717         9431 :          t402 = t108*br_a1
    1718         9431 :          t404 = t2*t5*t33
    1719         9431 :          t409 = t339*t83
    1720         9431 :          t420 = t106*br_a1
    1721         9431 :          t427 = t318/0.9e1_dp
    1722         9431 :          t428 = t395/0.3e1_dp
    1723         9431 :          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         9431 :                 + 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         9431 :                               *t233*t392/0.24e2_dp)*sx
    1731         9431 :          t443 = t6*t33
    1732         9431 :          t444 = t443*gamma
    1733         9431 :          t446 = t3*t444*t317
    1734         9431 :          t450 = t31*t42
    1735         9431 :          t451 = t450*gamma
    1736         9431 :          t454 = t40*t51
    1737         9431 :          t455 = t454*gamma
    1738         9431 :          t458 = t49*t61
    1739         9431 :          t459 = t458*gamma
    1740         9431 :          t462 = t59*t189
    1741         9431 :          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         9431 :                 /0.729e3_dp*t57*t463
    1745         9431 :          t468 = t23*t466*t83
    1746         9431 :          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         9431 :                 /0.729e3_dp*t79*t463
    1750         9431 :          t484 = t200*t102*R*t481
    1751         9431 :          t487 = t234*t235*t96
    1752         9431 :          t501 = gamma*t140
    1753         9431 :          t504 = 0.4e1_dp/0.9e1_dp*t3*t443*t501*t110
    1754         9431 :          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         9431 :                                                                  t468 + t506)*t95*t237
    1758         9431 :          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         9431 :                  /0.3e1_dp*t514)*t106
    1761         9431 :          t521 = t2*t6
    1762         9431 :          t525 = t143*t101
    1763         9431 :          t529 = t466*t83
    1764         9431 :          t540 = t420*t521
    1765         9431 :          t547 = 0.4e1_dp/0.9e1_dp*t446
    1766         9431 :          t548 = t514/0.3e1_dp
    1767         9431 :          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         9431 :                 *t116
    1774              :          e_tau = e_tau + (rho*t552*t123/0.8e1_dp - t119*t549*t123/0.8e1_dp - t310 &
    1775         9431 :                           *t233*t511/0.24e2_dp)*sx
    1776         9431 :          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         9431 :                 *t57*t462
    1780         9431 :          t580 = t23*t578*t83
    1781         9431 :          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         9431 :                 *t79*t462
    1785         9431 :          t596 = t200*t102*R*t593
    1786         9431 :          t612 = t3*t6
    1787         9431 :          t613 = t33*t140
    1788         9431 :          t614 = t613*t110
    1789         9431 :          t616 = t612*t614/0.9e1_dp
    1790         9431 :          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         9431 :                                                                  - t580 + t618)*t95*t237
    1794         9431 :          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         9431 :                  /0.3e1_dp*t626)*t106
    1797         9431 :          t638 = t578*t83
    1798         9431 :          t654 = t566/0.9e1_dp
    1799         9431 :          t655 = t626/0.3e1_dp
    1800         9431 :          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         9431 :                 + 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         9431 :                                           *t233*t623/0.24e2_dp)*sx
    1808              :       END IF
    1809              : 
    1810        17545 :    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     16977565 :    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     16977565 :       t1 = 8**(0.1e1_dp/0.3e1_dp)
    1844     16977565 :       t2 = t1**2
    1845     16977565 :       t3 = 0.1e1_dp/br_BB
    1846     16977565 :       t4 = pi**(0.1e1_dp/0.3e1_dp)
    1847     16977565 :       t5 = t4**2
    1848     16977565 :       t6 = 0.1e1_dp/t5
    1849     16977565 :       t8 = rho**(0.1e1_dp/0.3e1_dp)
    1850     16977565 :       t9 = t8**2
    1851     16977565 :       t10 = t9*rho
    1852     16977565 :       t11 = 0.1e1_dp/t10
    1853     16977565 :       t14 = ndrho**2
    1854     16977565 :       t15 = 0.1e1_dp/rho
    1855     16977565 :       t21 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t14*t15/0.4e1_dp)/0.3e1_dp
    1856     16977565 :       t25 = br_BB**2
    1857     16977565 :       t26 = t4*pi
    1858     16977565 :       t28 = rho**2
    1859     16977565 :       t31 = t21**2
    1860     16977565 :       t33 = t8*t28*rho/t31
    1861     16977565 :       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     16977565 :                 *t37*t3*t6*t11*t21)
    1864     16977565 :       t45 = t44 + 0.2e1_dp
    1865     16977565 :       t46 = t45**2
    1866     16977565 :       t50 = t10/t21
    1867     16977565 :       t56 = pi**2
    1868     16977565 :       t58 = t28**2
    1869     16977565 :       t62 = t58*rho/t31/t21
    1870     16977565 :       t65 = t5*t56
    1871     16977565 :       t69 = t31**2
    1872     16977565 :       t71 = t9*t58*t28/t69
    1873     16977565 :       t75 = t4*t56*pi
    1874     16977565 :       t77 = t58**2
    1875     16977565 :       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     16977565 :             *br_d4*t65*t71 + 0.32e2_dp/0.243e3_dp*br_d5*t75*t81
    1879     16977565 :       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     16977565 :              *br_e4*t65*t71 + 0.32e2_dp/0.243e3_dp*br_e5*t75*t81
    1883     16977565 :       t104 = t103**2
    1884     16977565 :       t111 = EXP(-t45*t84/t103)
    1885              :       t116 = (t46*t45*t85*t84/t104/t103*t111/0.3141592654e1_dp &
    1886     16977565 :               *t15)**(0.1e1_dp/0.3e1_dp)
    1887     16977565 :       brval = REAL(t2, KIND=dp)*t116/0.8e1_dp
    1888              : 
    1889     16977565 :       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       804729 :                                             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     16172836 :                                              sx, R, gamma, grad_deriv)
    1897              :       END IF
    1898              : 
    1899     16977565 :    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       804729 :    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       804729 :       IF (grad_deriv >= 0) THEN
    1947       804729 :          t1 = 0.1e1_dp/br_BB
    1948       804729 :          t2 = pi**(0.1e1_dp/0.3e1_dp)
    1949       804729 :          t3 = t2**2
    1950       804729 :          t4 = 0.1e1_dp/t3
    1951       804729 :          t5 = t1*t4
    1952       804729 :          t6 = rho**(0.1e1_dp/0.3e1_dp)
    1953       804729 :          t7 = t6**2
    1954       804729 :          t8 = t7*rho
    1955       804729 :          t9 = 0.1e1_dp/t8
    1956       804729 :          t12 = ndrho**2
    1957       804729 :          t13 = 0.1e1_dp/rho
    1958       804729 :          t19 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t12*t13/0.4e1_dp)/0.3e1_dp
    1959       804729 :          t20 = t9*t19
    1960       804729 :          t23 = br_BB**2
    1961       804729 :          t24 = t2*pi
    1962       804729 :          t25 = t23*t24
    1963       804729 :          t26 = rho**2
    1964       804729 :          t27 = t26*rho
    1965       804729 :          t28 = t6*t27
    1966       804729 :          t29 = t19**2
    1967       804729 :          t30 = 0.1e1_dp/t29
    1968       804729 :          t31 = t28*t30
    1969       804729 :          t35 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t31)
    1970       804729 :          t36 = t35*t1
    1971       804729 :          t37 = t4*t9
    1972              :          t41 = 0.1500000000000000e1_dp*t5*t20 + 0.3e1_dp/0.2e1_dp*t36*t37* &
    1973       804729 :                t19
    1974       804729 :          t42 = LOG(t41)
    1975       804729 :          t43 = t42 + 0.2e1_dp
    1976       804729 :          t44 = br_d1*t3
    1977       804729 :          t45 = 0.1e1_dp/t19
    1978       804729 :          t46 = t8*t45
    1979       804729 :          t49 = br_d2*t24
    1980       804729 :          t52 = pi**2
    1981       804729 :          t53 = br_d3*t52
    1982       804729 :          t54 = t26**2
    1983       804729 :          t55 = t54*rho
    1984       804729 :          t57 = 0.1e1_dp/t29/t19
    1985       804729 :          t58 = t55*t57
    1986       804729 :          t61 = t3*t52
    1987       804729 :          t62 = br_d4*t61
    1988       804729 :          t63 = t54*t26
    1989       804729 :          t64 = t7*t63
    1990       804729 :          t65 = t29**2
    1991       804729 :          t66 = 0.1e1_dp/t65
    1992       804729 :          t67 = t64*t66
    1993       804729 :          t71 = t2*t52*pi
    1994       804729 :          t72 = br_d5*t71
    1995       804729 :          t73 = t54**2
    1996       804729 :          t74 = t6*t73
    1997       804729 :          t76 = 0.1e1_dp/t65/t19
    1998       804729 :          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       804729 :                /0.243e3_dp*t72*t77
    2002       804729 :          t81 = t43*t80
    2003       804729 :          t82 = br_e1*t3
    2004       804729 :          t85 = br_e2*t24
    2005       804729 :          t88 = br_e3*t52
    2006       804729 :          t91 = br_e4*t61
    2007       804729 :          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       804729 :                /0.243e3_dp*t94*t77
    2011       804729 :          t98 = 0.1e1_dp/t97
    2012       804729 :          t99 = t81*t98
    2013       804729 :          t100 = 8**(0.1e1_dp/0.3e1_dp)
    2014       804729 :          t101 = t43**2
    2015       804729 :          t102 = t101*t43
    2016       804729 :          t103 = t80**2
    2017       804729 :          t104 = t103*t80
    2018       804729 :          t105 = t102*t104
    2019       804729 :          t106 = t97**2
    2020       804729 :          t108 = 0.1e1_dp/t106/t97
    2021       804729 :          t109 = t105*t108
    2022       804729 :          t110 = EXP(-t99)
    2023       804729 :          t111 = 0.1e1_dp/0.3141592654e1_dp
    2024       804729 :          t112 = t110*t111
    2025       804729 :          t114 = t109*t112*t13
    2026       804729 :          t115 = t114**(0.1e1_dp/0.3e1_dp)
    2027       804729 :          t116 = 0.1e1_dp/t115
    2028       804729 :          t117 = REAL(t100, KIND=dp)*t116
    2029       804729 :          t118 = t117*R
    2030       804729 :          t119 = t99*t118
    2031       804729 :          t121 = EXP(t99 - t119)
    2032       804729 :          t123 = t121*t43
    2033       804729 :          t124 = t80*t98
    2034       804729 :          t125 = t123*t124
    2035       804729 :          t129 = t98*REAL(t100, KIND=dp)*t116*R
    2036       804729 :          t132 = EXP(-t99 - t119)
    2037       804729 :          t134 = t132*t43
    2038       804729 :          t135 = t134*t124
    2039              :          t138 = -0.2e1_dp*t121 + t125 - t123*t80*t129 + 0.2e1_dp*t132 + t135 &
    2040       804729 :                 + t134*t80*t129
    2041       804729 :          t139 = rho*t138
    2042       804729 :          e_0 = e_0 + (t139*t117/0.8e1_dp)*sx
    2043              :       END IF
    2044       804729 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    2045       415000 :          t145 = 0.1e1_dp/t7/t26
    2046       415000 :          t152 = 0.1e1_dp/t7/t27*gamma*t12
    2047       415000 :          t155 = 0.1e1_dp/t35
    2048       415000 :          t158 = t6*t26
    2049       415000 :          t159 = t158*t30
    2050       415000 :          t162 = t6*rho
    2051       415000 :          t164 = t57*gamma
    2052       415000 :          t165 = t164*t12
    2053       415000 :          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       415000 :                 *t4*t145*t19 - t176*t152/0.8e1_dp
    2058       415000 :          t180 = 0.1e1_dp/t41
    2059       415000 :          t181 = t179*t180
    2060       415000 :          t182 = t181*t124
    2061       415000 :          t183 = t7*t45
    2062       415000 :          t186 = 0.1e1_dp/t6
    2063       415000 :          t188 = t30*gamma
    2064       415000 :          t189 = t188*t12
    2065       415000 :          t197 = t54*t57
    2066       415000 :          t201 = t66*gamma
    2067       415000 :          t202 = t201*t12
    2068       415000 :          t205 = t7*t55
    2069       415000 :          t206 = t205*t66
    2070       415000 :          t209 = t7*t54
    2071       415000 :          t211 = t76*gamma
    2072       415000 :          t212 = t211*t12
    2073       415000 :          t216 = t6*t54*t27
    2074       415000 :          t217 = t216*t76
    2075       415000 :          t220 = t6*t63
    2076       415000 :          t223 = 0.1e1_dp/t65/t29
    2077       415000 :          t224 = t223*gamma
    2078       415000 :          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       415000 :                 t220*t225
    2085       415000 :          t230 = t43*t228*t98
    2086       415000 :          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       415000 :                 t220*t225
    2093       415000 :          t259 = t81*t231*t257
    2094       415000 :          t261 = t181*t80*t129
    2095       415000 :          t262 = t230*t118
    2096       415000 :          t263 = t81*t231
    2097       415000 :          t265 = t117*R*t257
    2098       415000 :          t266 = t263*t265
    2099       415000 :          t269 = REAL(t100, KIND=dp)/t115/t114
    2100       415000 :          t272 = t101*t104*t108*t110
    2101       415000 :          t273 = t111*t13
    2102       415000 :          t278 = t102*t103*t108
    2103       415000 :          t283 = t106**2
    2104       415000 :          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       415000 :                 *t110*t273 - t109*t112/t26
    2108       415000 :          t299 = t269*R*t297
    2109       415000 :          t301 = t99*t299/0.3e1_dp
    2110       415000 :          t303 = (t182 + t230 - t259 - t261 - t262 + t266 + t301)*t121
    2111       415000 :          t306 = t121*t179
    2112       415000 :          t307 = t180*t80
    2113       415000 :          t308 = t307*t98
    2114       415000 :          t310 = t228*t98
    2115       415000 :          t312 = t80*t231
    2116       415000 :          t313 = t312*t257
    2117       415000 :          t321 = t123*t312
    2118       415000 :          t326 = (-t182 - t230 + t259 - t261 - t262 + t266 + t301)*t132
    2119       415000 :          t329 = t132*t179
    2120       415000 :          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       415000 :                 /0.3e1_dp
    2127              :          e_rho = e_rho + (t138*REAL(t100, KIND=dp)*t116/0.8e1_dp + rho*t343*t117/0.8e1_dp &
    2128       415000 :                           - t139*t269*t297/0.24e2_dp)*sx
    2129       415000 :          t351 = t145*gamma*ndrho
    2130       415000 :          t354 = t155*br_BB
    2131       415000 :          t355 = t354*t3
    2132              :          t363 = 0.2500000000000000e0_dp*t5*t351 - t355*t7*t30*gamma*ndrho &
    2133       415000 :                 /0.9e1_dp + t176*t351/0.4e1_dp
    2134       415000 :          t364 = t363*t180
    2135       415000 :          t365 = t364*t124
    2136       415000 :          t367 = t188*ndrho
    2137       415000 :          t371 = t164*ndrho
    2138       415000 :          t375 = t201*ndrho
    2139       415000 :          t379 = t211*ndrho
    2140       415000 :          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       415000 :                 *t379 - 0.80e2_dp/0.729e3_dp*t72*t216*t383
    2144       415000 :          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       415000 :                 *t379 - 0.80e2_dp/0.729e3_dp*t94*t216*t383
    2148       415000 :          t406 = t81*t231*t404
    2149       415000 :          t408 = t364*t80*t129
    2150       415000 :          t409 = t388*t118
    2151       415000 :          t411 = t117*R*t404
    2152       415000 :          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       415000 :                 *t110*t273
    2156       415000 :          t430 = t269*R*t428
    2157       415000 :          t432 = t99*t430/0.3e1_dp
    2158       415000 :          t434 = (t365 + t388 - t406 - t408 - t409 + t412 + t432)*t121
    2159       415000 :          t437 = t121*t363
    2160       415000 :          t439 = t386*t98
    2161       415000 :          t441 = t312*t404
    2162       415000 :          t453 = (-t365 - t388 + t406 - t408 - t409 + t412 + t432)*t132
    2163       415000 :          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       415000 :                 /0.3e1_dp
    2170       415000 :          e_ndrho = e_ndrho + (rho*t469*t117/0.8e1_dp - t139*t269*t428/0.24e2_dp)*sx
    2171       415000 :          t479 = t8*t30
    2172       415000 :          t480 = t479*gamma
    2173              :          t485 = -0.1000000000e1_dp*t5*t9*gamma + 0.4e1_dp/0.9e1_dp*t355*t480 &
    2174       415000 :                 - t36*t37*gamma
    2175       415000 :          t486 = t485*t180
    2176       415000 :          t487 = t486*t124
    2177       415000 :          t490 = t28*t57
    2178       415000 :          t491 = t490*gamma
    2179       415000 :          t494 = t55*t66
    2180       415000 :          t495 = t494*gamma
    2181       415000 :          t498 = t64*t76
    2182       415000 :          t499 = t498*gamma
    2183       415000 :          t502 = t74*t223
    2184       415000 :          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       415000 :                 /0.729e3_dp*t72*t503
    2188       415000 :          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       415000 :                 /0.729e3_dp*t94*t503
    2192       415000 :          t521 = t81*t231*t519
    2193       415000 :          t523 = t486*t80*t129
    2194       415000 :          t524 = t508*t118
    2195       415000 :          t526 = t117*R*t519
    2196       415000 :          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       415000 :                 *t110*t273
    2200       415000 :          t545 = t269*R*t543
    2201       415000 :          t547 = t99*t545/0.3e1_dp
    2202       415000 :          t549 = (t487 + t508 - t521 - t523 - t524 + t527 + t547)*t121
    2203       415000 :          t552 = t121*t485
    2204       415000 :          t554 = t506*t98
    2205       415000 :          t556 = t312*t519
    2206       415000 :          t568 = (-t487 - t508 + t521 - t523 - t524 + t527 + t547)*t132
    2207       415000 :          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       415000 :                 /0.3e1_dp
    2214       415000 :          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       415000 :                 + t36*t37/0.4e1_dp
    2217       415000 :          t600 = t599*t180
    2218       415000 :          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       415000 :                 *t72*t502
    2222       415000 :          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       415000 :                 *t94*t502
    2226       415000 :          t627 = t81*t231*t625
    2227       415000 :          t629 = t600*t80*t129
    2228       415000 :          t630 = t614*t118
    2229       415000 :          t632 = t117*R*t625
    2230       415000 :          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       415000 :                 *t110*t273
    2234       415000 :          t651 = t269*R*t649
    2235       415000 :          t653 = t99*t651/0.3e1_dp
    2236       415000 :          t655 = (t601 + t614 - t627 - t629 - t630 + t633 + t653)*t121
    2237       415000 :          t658 = t121*t599
    2238       415000 :          t660 = t612*t98
    2239       415000 :          t662 = t312*t625
    2240       415000 :          t674 = (-t601 - t614 + t627 - t629 - t630 + t633 + t653)*t132
    2241       415000 :          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       415000 :                 /0.3e1_dp
    2248       415000 :          e_laplace_rho = e_laplace_rho + (rho*t690*t117/0.8e1_dp - t139*t269*t649/0.24e2_dp)*sx
    2249              :       END IF
    2250              : 
    2251       804729 :    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     16172836 :    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     16172836 :       IF (grad_deriv >= 0) THEN
    2298     16172836 :          t1 = 0.1e1_dp/br_BB
    2299     16172836 :          t2 = pi**(0.1e1_dp/0.3e1_dp)
    2300     16172836 :          t3 = t2**2
    2301     16172836 :          t4 = 0.1e1_dp/t3
    2302     16172836 :          t5 = t1*t4
    2303     16172836 :          t6 = rho**(0.1e1_dp/0.3e1_dp)
    2304     16172836 :          t7 = t6**2
    2305     16172836 :          t8 = t7*rho
    2306     16172836 :          t9 = 0.1e1_dp/t8
    2307     16172836 :          t12 = ndrho**2
    2308     16172836 :          t13 = 0.1e1_dp/rho
    2309     16172836 :          t19 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t12*t13/0.4e1_dp)/0.3e1_dp
    2310     16172836 :          t20 = t9*t19
    2311     16172836 :          t23 = br_BB**2
    2312     16172836 :          t24 = t2*pi
    2313     16172836 :          t25 = t23*t24
    2314     16172836 :          t26 = rho**2
    2315     16172836 :          t27 = t26*rho
    2316     16172836 :          t28 = t6*t27
    2317     16172836 :          t29 = t19**2
    2318     16172836 :          t30 = 0.1e1_dp/t29
    2319     16172836 :          t31 = t28*t30
    2320     16172836 :          t35 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t31)
    2321     16172836 :          t36 = t35*t1
    2322     16172836 :          t37 = t4*t9
    2323              :          t41 = 0.1500000000000000e1_dp*t5*t20 + 0.3e1_dp/0.2e1_dp*t36*t37* &
    2324     16172836 :                t19
    2325     16172836 :          t42 = LOG(t41)
    2326     16172836 :          t43 = t42 + 0.2e1_dp
    2327     16172836 :          t44 = br_d1*t3
    2328     16172836 :          t45 = 0.1e1_dp/t19
    2329     16172836 :          t46 = t8*t45
    2330     16172836 :          t49 = br_d2*t24
    2331     16172836 :          t52 = pi**2
    2332     16172836 :          t53 = br_d3*t52
    2333     16172836 :          t54 = t26**2
    2334     16172836 :          t55 = t54*rho
    2335     16172836 :          t57 = 0.1e1_dp/t29/t19
    2336     16172836 :          t58 = t55*t57
    2337     16172836 :          t61 = t3*t52
    2338     16172836 :          t62 = br_d4*t61
    2339     16172836 :          t63 = t54*t26
    2340     16172836 :          t64 = t7*t63
    2341     16172836 :          t65 = t29**2
    2342     16172836 :          t66 = 0.1e1_dp/t65
    2343     16172836 :          t67 = t64*t66
    2344     16172836 :          t71 = t2*t52*pi
    2345     16172836 :          t72 = br_d5*t71
    2346     16172836 :          t73 = t54**2
    2347     16172836 :          t74 = t6*t73
    2348     16172836 :          t76 = 0.1e1_dp/t65/t19
    2349     16172836 :          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     16172836 :                /0.243e3_dp*t72*t77
    2353     16172836 :          t81 = t43*t80
    2354     16172836 :          t82 = br_e1*t3
    2355     16172836 :          t85 = br_e2*t24
    2356     16172836 :          t88 = br_e3*t52
    2357     16172836 :          t91 = br_e4*t61
    2358     16172836 :          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     16172836 :                /0.243e3_dp*t94*t77
    2362     16172836 :          t98 = 0.1e1_dp/t97
    2363     16172836 :          t99 = t81*t98
    2364     16172836 :          t100 = 8**(0.1e1_dp/0.3e1_dp)
    2365     16172836 :          t101 = t43**2
    2366     16172836 :          t102 = t101*t43
    2367     16172836 :          t103 = t80**2
    2368     16172836 :          t104 = t103*t80
    2369     16172836 :          t105 = t102*t104
    2370     16172836 :          t106 = t97**2
    2371     16172836 :          t108 = 0.1e1_dp/t106/t97
    2372     16172836 :          t109 = t105*t108
    2373     16172836 :          t110 = EXP(-t99)
    2374     16172836 :          t111 = 0.1e1_dp/0.3141592654e1_dp
    2375     16172836 :          t112 = t110*t111
    2376     16172836 :          t114 = t109*t112*t13
    2377     16172836 :          t115 = t114**(0.1e1_dp/0.3e1_dp)
    2378     16172836 :          t116 = 0.1e1_dp/t115
    2379     16172836 :          t117 = REAL(t100, KIND=dp)*t116
    2380     16172836 :          t118 = t117*R
    2381     16172836 :          t119 = t99*t118
    2382     16172836 :          t121 = EXP(0.2e1_dp*t119)
    2383     16172836 :          t123 = t121*R
    2384     16172836 :          t124 = t123*t43
    2385     16172836 :          t125 = t80*t98
    2386     16172836 :          t126 = t125*t117
    2387     16172836 :          t128 = t121*t43
    2388     16172836 :          t130 = t99 + t119
    2389     16172836 :          t131 = EXP(t130)
    2390              :          t133 = 0.2e1_dp*t121 - t124*t126 + t128*t125 + 0.2e1_dp + t99 + t119 &
    2391     16172836 :                 - 0.4e1_dp*t131
    2392     16172836 :          t134 = rho*t133
    2393     16172836 :          t136 = EXP(-t130)
    2394     16172836 :          t138 = t136*REAL(t100, KIND=dp)*t116
    2395     16172836 :          e_0 = e_0 + (t134*t138/0.8e1_dp)*sx
    2396              :       END IF
    2397     16172836 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    2398      8506079 :          t144 = 0.1e1_dp/t7/t26
    2399      8506079 :          t151 = 0.1e1_dp/t7/t27*gamma*t12
    2400      8506079 :          t154 = 0.1e1_dp/t35
    2401      8506079 :          t157 = t6*t26
    2402      8506079 :          t158 = t157*t30
    2403      8506079 :          t161 = t6*rho
    2404      8506079 :          t163 = t57*gamma
    2405      8506079 :          t164 = t163*t12
    2406      8506079 :          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      8506079 :                 *t4*t144*t19 - t175*t151/0.8e1_dp
    2411      8506079 :          t179 = 0.1e1_dp/t41
    2412      8506079 :          t180 = t178*t179
    2413      8506079 :          t182 = t98*REAL(t100, KIND=dp)
    2414      8506079 :          t184 = t182*t116*R
    2415      8506079 :          t185 = t180*t80*t184
    2416      8506079 :          t187 = t7*t45
    2417      8506079 :          t190 = 0.1e1_dp/t6
    2418      8506079 :          t192 = t30*gamma
    2419      8506079 :          t193 = t192*t12
    2420      8506079 :          t201 = t54*t57
    2421      8506079 :          t205 = t66*gamma
    2422      8506079 :          t206 = t205*t12
    2423      8506079 :          t209 = t7*t55
    2424      8506079 :          t210 = t209*t66
    2425      8506079 :          t213 = t7*t54
    2426      8506079 :          t215 = t76*gamma
    2427      8506079 :          t216 = t215*t12
    2428      8506079 :          t220 = t6*t54*t27
    2429      8506079 :          t221 = t220*t76
    2430      8506079 :          t224 = t6*t63
    2431      8506079 :          t227 = 0.1e1_dp/t65/t29
    2432      8506079 :          t228 = t227*gamma
    2433      8506079 :          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      8506079 :                 t224*t229
    2440      8506079 :          t234 = t43*t232*t98
    2441      8506079 :          t235 = t234*t118
    2442      8506079 :          t237 = 0.1e1_dp/t106
    2443      8506079 :          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      8506079 :                 t224*t229
    2450      8506079 :          t267 = t238*t117*R*t264
    2451      8506079 :          t270 = 0.1e1_dp/t115/t114
    2452      8506079 :          t271 = REAL(t100, KIND=dp)*t270
    2453      8506079 :          t274 = t101*t104*t108*t110
    2454      8506079 :          t275 = t111*t13
    2455      8506079 :          t280 = t102*t103*t108
    2456      8506079 :          t285 = t106**2
    2457      8506079 :          t287 = t105/t285
    2458      8506079 :          t292 = t180*t125
    2459      8506079 :          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      8506079 :                 *t110*t275 - t109*t112/t26
    2463      8506079 :          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      8506079 :                  *t305)*t121
    2466      8506079 :          t310 = R*t43
    2467      8506079 :          t315 = t232*t98
    2468      8506079 :          t318 = t123*t81
    2469      8506079 :          t319 = t237*REAL(t100, KIND=dp)
    2470      8506079 :          t330 = t179*t80*t98
    2471      8506079 :          t333 = t80*t237
    2472      8506079 :          t336 = t305/0.3e1_dp
    2473      8506079 :          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      8506079 :                 - 0.4e1_dp*t337*t131
    2479      8506079 :          t348 = t134*t136
    2480              :          e_rho = e_rho + (t133*t136*t117/0.8e1_dp + rho*t340*t138/0.8e1_dp - t134 &
    2481      8506079 :                           *t337*t138/0.8e1_dp - t348*t271*t302/0.24e2_dp)*sx
    2482      8506079 :          t353 = t144*gamma*ndrho
    2483      8506079 :          t356 = t154*br_BB
    2484      8506079 :          t357 = t356*t3
    2485              :          t365 = 0.2500000000000000e0_dp*t5*t353 - t357*t7*t30*gamma*ndrho &
    2486      8506079 :                 /0.9e1_dp + t175*t353/0.4e1_dp
    2487      8506079 :          t366 = t365*t179
    2488      8506079 :          t368 = t366*t80*t184
    2489      8506079 :          t371 = t192*ndrho
    2490      8506079 :          t375 = t163*ndrho
    2491      8506079 :          t379 = t205*ndrho
    2492      8506079 :          t383 = t215*ndrho
    2493      8506079 :          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      8506079 :                 *t383 - 0.80e2_dp/0.729e3_dp*t72*t220*t387
    2497      8506079 :          t392 = t43*t390*t98
    2498      8506079 :          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      8506079 :                 *t383 - 0.80e2_dp/0.729e3_dp*t94*t220*t387
    2502      8506079 :          t413 = t238*t117*R*t410
    2503      8506079 :          t426 = t366*t125
    2504      8506079 :          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      8506079 :                 *t110*t275
    2508      8506079 :          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      8506079 :                  *t436)*t121
    2511      8506079 :          t445 = t390*t98
    2512      8506079 :          t461 = t436/0.3e1_dp
    2513      8506079 :          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      8506079 :                 - 0.4e1_dp*t462*t131
    2519              :          e_ndrho = e_ndrho + (rho*t465*t138/0.8e1_dp - t134*t462*t138/0.8e1_dp - t348 &
    2520      8506079 :                               *t271*t433/0.24e2_dp)*sx
    2521      8506079 :          t479 = t8*t30
    2522      8506079 :          t480 = t479*gamma
    2523              :          t485 = -0.1000000000e1_dp*t5*t9*gamma + 0.4e1_dp/0.9e1_dp*t357*t480 &
    2524      8506079 :                 - t36*t37*gamma
    2525      8506079 :          t486 = t485*t179
    2526      8506079 :          t488 = t486*t80*t184
    2527      8506079 :          t492 = t28*t57
    2528      8506079 :          t493 = t492*gamma
    2529      8506079 :          t496 = t55*t66
    2530      8506079 :          t497 = t496*gamma
    2531      8506079 :          t500 = t64*t76
    2532      8506079 :          t501 = t500*gamma
    2533      8506079 :          t504 = t74*t227
    2534      8506079 :          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      8506079 :                 /0.729e3_dp*t72*t505
    2538      8506079 :          t510 = t43*t508*t98
    2539      8506079 :          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      8506079 :                 /0.729e3_dp*t94*t505
    2543      8506079 :          t526 = t238*t117*R*t523
    2544      8506079 :          t539 = t486*t125
    2545      8506079 :          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      8506079 :                 *t110*t275
    2549      8506079 :          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      8506079 :                  *t549)*t121
    2552      8506079 :          t558 = t508*t98
    2553      8506079 :          t574 = t549/0.3e1_dp
    2554      8506079 :          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      8506079 :                 - 0.4e1_dp*t575*t131
    2560              :          e_tau = e_tau + (rho*t578*t138/0.8e1_dp - t134*t575*t138/0.8e1_dp - t348 &
    2561      8506079 :                           *t271*t546/0.24e2_dp)*sx
    2562              :          t597 = 0.2500000000000000e0_dp*t5*t9 - t356*t3*t8*t30/0.9e1_dp &
    2563      8506079 :                 + t36*t37/0.4e1_dp
    2564      8506079 :          t598 = t597*t179
    2565      8506079 :          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      8506079 :                 *t72*t504
    2569      8506079 :          t614 = t43*t612*t98
    2570      8506079 :          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      8506079 :                 *t94*t504
    2574      8506079 :          t630 = t238*t117*R*t627
    2575      8506079 :          t643 = t598*t125
    2576      8506079 :          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      8506079 :                 *t110*t275
    2580      8506079 :          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      8506079 :                  *t653)*t121
    2583      8506079 :          t662 = t612*t98
    2584      8506079 :          t678 = t653/0.3e1_dp
    2585      8506079 :          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      8506079 :                 - 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      8506079 :                                           *t271*t650/0.24e2_dp)*sx
    2593              :       END IF
    2594              : 
    2595     16172836 :    END SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_lte_b
    2596              : 
    2597              : END MODULE xc_xbecke_roussel
        

Generated by: LCOV version 2.0-1