LCOV - code coverage report
Current view: top level - src/xc - xc_xbecke88_long_range.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 61.2 % 771 472
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 6 6

            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 longrange part of Becke 88 exchange functional
      10              : !> \par History
      11              : !>      01.2008 created [mguidon]
      12              : !> \author Manuel Guidon
      13              : ! **************************************************************************************************
      14              : MODULE xc_xbecke88_long_range
      15              :    USE bibliography,                    ONLY: Becke1988,&
      16              :                                               cite_reference
      17              :    USE cp_array_utils,                  ONLY: cp_3d_r_cp_type
      18              :    USE input_section_types,             ONLY: section_vals_type,&
      19              :                                               section_vals_val_get
      20              :    USE kinds,                           ONLY: dp
      21              :    USE mathconstants,                   ONLY: pi,&
      22              :                                               rootpi
      23              :    USE xc_derivative_desc,              ONLY: deriv_norm_drho,&
      24              :                                               deriv_norm_drhoa,&
      25              :                                               deriv_norm_drhob,&
      26              :                                               deriv_rho,&
      27              :                                               deriv_rhoa,&
      28              :                                               deriv_rhob
      29              :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      30              :                                               xc_dset_get_derivative
      31              :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      32              :                                               xc_derivative_type
      33              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      34              :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      35              :                                               xc_rho_set_type
      36              : #include "../base/base_uses.f90"
      37              : 
      38              :    IMPLICIT NONE
      39              :    PRIVATE
      40              : 
      41              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      42              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbecke88_long_range'
      43              :    REAL(kind=dp), PARAMETER :: beta = 0.0042_dp
      44              : 
      45              :    PUBLIC :: xb88_lr_lda_info, xb88_lr_lsd_info, xb88_lr_lda_eval, xb88_lr_lsd_eval
      46              : CONTAINS
      47              : 
      48              : ! **************************************************************************************************
      49              : !> \brief return various information on the functional
      50              : !> \param reference string with the reference of the actual functional
      51              : !> \param shortform string with the shortform of the functional name
      52              : !> \param needs the components needed by this functional are set to
      53              : !>        true (does not set the unneeded components to false)
      54              : !> \param max_deriv ...
      55              : !> \par History
      56              : !>      01.2008 created [mguidon]
      57              : !> \author Manuel Guidon
      58              : ! **************************************************************************************************
      59         1316 :    SUBROUTINE xb88_lr_lda_info(reference, shortform, needs, max_deriv)
      60              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      61              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      62              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      63              : 
      64         1316 :       IF (PRESENT(reference)) THEN
      65            8 :          reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LDA version}"
      66              :       END IF
      67         1316 :       IF (PRESENT(shortform)) THEN
      68            8 :          shortform = "Becke 1988 Exchange Functional (LDA)"
      69              :       END IF
      70         1316 :       IF (PRESENT(needs)) THEN
      71         1308 :          needs%rho = .TRUE.
      72         1308 :          needs%norm_drho = .TRUE.
      73              :       END IF
      74         1316 :       IF (PRESENT(max_deriv)) max_deriv = 3
      75              : 
      76         1316 :    END SUBROUTINE xb88_lr_lda_info
      77              : 
      78              : ! **************************************************************************************************
      79              : !> \brief return various information on the functional
      80              : !> \param reference string with the reference of the actual functional
      81              : !> \param shortform string with the shortform of the functional name
      82              : !> \param needs the components needed by this functional are set to
      83              : !>        true (does not set the unneeded components to false)
      84              : !> \param max_deriv ...
      85              : !> \par History
      86              : !>      01.2008 created [mguidon]
      87              : !> \author Manuel Guidon
      88              : ! **************************************************************************************************
      89           61 :    SUBROUTINE xb88_lr_lsd_info(reference, shortform, needs, max_deriv)
      90              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      91              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      92              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      93              : 
      94           61 :       IF (PRESENT(reference)) THEN
      95            1 :          reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LSD version}"
      96              :       END IF
      97           61 :       IF (PRESENT(shortform)) THEN
      98            1 :          shortform = "Becke 1988 Exchange Functional (LSD)"
      99              :       END IF
     100           61 :       IF (PRESENT(needs)) THEN
     101           60 :          needs%rho_spin = .TRUE.
     102           60 :          needs%rho_spin_1_3 = .TRUE.
     103           60 :          needs%norm_drho_spin = .TRUE.
     104              :       END IF
     105           61 :       IF (PRESENT(max_deriv)) max_deriv = 3
     106              : 
     107           61 :    END SUBROUTINE xb88_lr_lsd_info
     108              : 
     109              : ! **************************************************************************************************
     110              : !> \brief evaluates the becke 88 longrange exchange functional for lda
     111              : !> \param rho_set the density where you want to evaluate the functional
     112              : !> \param deriv_set place where to store the functional derivatives (they are
     113              : !>        added to the derivatives)
     114              : !> \param grad_deriv degree of the derivative that should be evaluated,
     115              : !>        if positive all the derivatives up to the given degree are evaluated,
     116              : !>        if negative only the given degree is calculated
     117              : !> \param xb88_lr_params input parameters (scaling)
     118              : !> \par History
     119              : !>      01.2008 created [mguidon]
     120              : !> \author Manuel Guidon
     121              : ! **************************************************************************************************
     122         4888 :    SUBROUTINE xb88_lr_lda_eval(rho_set, deriv_set, grad_deriv, xb88_lr_params)
     123              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     124              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     125              :       INTEGER, INTENT(in)                                :: grad_deriv
     126              :       TYPE(section_vals_type), POINTER                   :: xb88_lr_params
     127              : 
     128              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xb88_lr_lda_eval'
     129              : 
     130              :       INTEGER                                            :: handle, npoints
     131              :       INTEGER, DIMENSION(2, 3)                           :: bo
     132              :       REAL(kind=dp)                                      :: epsilon_rho, omega, sx
     133         1222 :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
     134         1222 :          e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
     135         1222 :          e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
     136              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     137              : 
     138         1222 :       CALL timeset(routineN, handle)
     139              : 
     140         1222 :       CALL section_vals_val_get(xb88_lr_params, "scale_x", r_val=sx)
     141         1222 :       CALL section_vals_val_get(xb88_lr_params, "omega", r_val=omega)
     142              : 
     143         1222 :       CALL cite_reference(Becke1988)
     144              : 
     145              :       CALL xc_rho_set_get(rho_set, rho=rho, &
     146         1222 :                           norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
     147         1222 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     148              : 
     149              :       ! meaningful default for the arrays we don't need: let us make compiler
     150              :       ! and debugger happy...
     151         1222 :       dummy => rho
     152              : 
     153         1222 :       e_0 => dummy
     154         1222 :       e_rho => dummy
     155         1222 :       e_ndrho => dummy
     156         1222 :       e_rho_rho => dummy
     157         1222 :       e_ndrho_rho => dummy
     158         1222 :       e_ndrho_ndrho => dummy
     159         1222 :       e_rho_rho_rho => dummy
     160         1222 :       e_ndrho_rho_rho => dummy
     161         1222 :       e_ndrho_ndrho_rho => dummy
     162         1222 :       e_ndrho_ndrho_ndrho => dummy
     163              : 
     164         1222 :       IF (grad_deriv >= 0) THEN
     165              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     166         1222 :                                          allocate_deriv=.TRUE.)
     167         1222 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     168              :       END IF
     169         1222 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     170              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     171         1208 :                                          allocate_deriv=.TRUE.)
     172         1208 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     173              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     174         1208 :                                          allocate_deriv=.TRUE.)
     175         1208 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     176              :       END IF
     177         1222 :       IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     178              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     179          348 :                                          allocate_deriv=.TRUE.)
     180          348 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     181              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
     182          348 :                                          allocate_deriv=.TRUE.)
     183          348 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
     184              :          deriv => xc_dset_get_derivative(deriv_set, &
     185          348 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     186          348 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     187              :       END IF
     188         1222 :       IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
     189              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
     190            0 :                                          allocate_deriv=.TRUE.)
     191            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     192              :          deriv => xc_dset_get_derivative(deriv_set, &
     193            0 :                                          [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.TRUE.)
     194            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
     195              :          deriv => xc_dset_get_derivative(deriv_set, &
     196            0 :                                          [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.TRUE.)
     197            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
     198              :          deriv => xc_dset_get_derivative(deriv_set, &
     199            0 :                                          [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     200            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
     201              :       END IF
     202         1222 :       IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
     203            0 :          CPABORT("derivatives bigger than 3 not implemented")
     204              :       END IF
     205              : 
     206              : !$OMP     PARALLEL DEFAULT(NONE) &
     207              : !$OMP              SHARED(rho, norm_drho, e_0, e_rho, e_ndrho, e_rho_rho) &
     208              : !$OMP              SHARED(e_ndrho_rho, e_ndrho_ndrho, e_rho_rho_rho) &
     209              : !$OMP              SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
     210              : !$OMP              SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
     211         1222 : !$OMP              SHARED(epsilon_rho, sx, omega)
     212              : 
     213              :       CALL xb88_lr_lda_calc(rho=rho, norm_drho=norm_drho, &
     214              :                             e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
     215              :                             e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
     216              :                             e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
     217              :                             e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
     218              :                             e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, grad_deriv=grad_deriv, &
     219              :                             npoints=npoints, epsilon_rho=epsilon_rho, &
     220              :                             sx=sx, omega=omega)
     221              : 
     222              : !$OMP     END PARALLEL
     223              : 
     224         1222 :       CALL timestop(handle)
     225         1222 :    END SUBROUTINE xb88_lr_lda_eval
     226              : 
     227              : ! **************************************************************************************************
     228              : !> \brief low level calculation of the becke 88 exchange functional for lda
     229              : !> \param rho alpha or beta spin density
     230              : !> \param norm_drho || grad rho ||
     231              : !> \param e_0 adds to it the local value of the functional
     232              : !> \param e_rho derivative of the functional wrt. to the variables
     233              : !>        named where the * is.
     234              : !> \param e_ndrho ...
     235              : !> \param e_rho_rho ...
     236              : !> \param e_ndrho_rho ...
     237              : !> \param e_ndrho_ndrho ...
     238              : !> \param e_rho_rho_rho ...
     239              : !> \param e_ndrho_rho_rho ...
     240              : !> \param e_ndrho_ndrho_rho ...
     241              : !> \param e_ndrho_ndrho_ndrho ...
     242              : !> \param grad_deriv ...
     243              : !> \param npoints ...
     244              : !> \param epsilon_rho ...
     245              : !> \param sx scaling-parameter for exchange
     246              : !> \param omega parameter for erfc
     247              : !> \par History
     248              : !>      01.2008 created [mguidon]
     249              : !> \author Manuel Guidon
     250              : !> \note
     251              : !>      - Just took the lsd code and scaled rho and ndrho by 1/2 (e_0 with 2.0)
     252              : !>      - Derivatives higher than 1 not tested
     253              : ! **************************************************************************************************
     254         1222 :    SUBROUTINE xb88_lr_lda_calc(rho, norm_drho, &
     255         1222 :                                e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
     256         1222 :                                e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
     257         1222 :                                e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, sx, &
     258              :                                omega)
     259              :       INTEGER, INTENT(in)                                :: npoints, grad_deriv
     260              :       REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_ndrho, &
     261              :          e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, &
     262              :          e_ndrho, e_rho, e_0
     263              :       REAL(kind=dp), DIMENSION(1:npoints), INTENT(in)    :: norm_drho, rho
     264              :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, sx, omega
     265              : 
     266              :       INTEGER                                            :: ii
     267              :       REAL(kind=dp) :: Cx, epsilon_rho43, my_epsilon_rho, my_ndrho, my_rho, t1, t10, t100, t1002, &
     268              :          t1009, t101, t1013, t103, t104, t1044, t105, t1069, t109, t1091, t11, t1102, t112, t113, &
     269              :          t1136, t1141, t1146, t1153, t1156, t116, t1163, t1167, t117, t1177, t118, t12, t122, &
     270              :          t1220, t126, t127, t128, t1284, t130, t132, t133, t1334, t1341, t1345, t135, t136, t137, &
     271              :          t1370, t1396, t140, t1400, t1404, t1405, t141, t1412, t143, t1449, t1456, t146, t1467, &
     272              :          t147, t1472, t148, t151, t155, t1553, t16, t168, t169, t17, t172, t173, t176, t177, t18, &
     273              :          t185, t186, t190, t196, t2, t200, t207, t21, t211, t212, t213
     274              :       REAL(kind=dp) :: t216, t219, t22, t221, t222, t225, t226, t23, t230, t231, t232, t233, t237, &
     275              :          t24, t241, t244, t245, t250, t251, t254, t255, t258, t259, t26, t264, t265, t27, t270, &
     276              :          t271, t28, t281, t284, t285, t289, t29, t293, t297, t3, t30, t304, t31, t311, t312, t313, &
     277              :          t316, t319, t321, t323, t325, t326, t328, t330, t335, t339, t34, t343, t346, t347, t351, &
     278              :          t354, t358, t36, t365, t368, t37, t372, t377, t38, t382, t383, t384, t39, t393, t397, &
     279              :          t40, t400, t401, t404, t405, t408, t41, t412, t413, t417, t418, t42, t428, t429, t43, &
     280              :          t435, t44, t451, t452, t455, t457, t46, t460, t462, t463, t464
     281              :       REAL(kind=dp) :: t465, t466, t467, t47, t470, t473, t478, t479, t48, t482, t486, t489, t49, &
     282              :          t492, t496, t5, t501, t502, t505, t508, t51, t513, t514, t519, t52, t521, t525, t526, &
     283              :          t529, t530, t533, t534, t536, t537, t539, t55, t562, t566, t570, t571, t572, t573, t574, &
     284              :          t575, t576, t580, t586, t59, t590, t6, t60, t601, t605, t606, t613, t624, t628, t632, &
     285              :          t639, t64, t640, t641, t657, t667, t669, t677, t68, t687, t69, t7, t71, t716, t72, t722, &
     286              :          t735, t739, t746, t75, t76, t769, t77, t79, t791, t793, t8, t838, t84, t842, t846, t85, &
     287              :          t857, t86, t860, t867, t87, t880, t90, t91, t933, t94, t95, t961
     288              :       REAL(kind=dp) :: t98, t99, xx
     289              : 
     290         1222 :       my_epsilon_rho = epsilon_rho
     291         1222 :       epsilon_rho43 = my_epsilon_rho**(4.0_dp/3.0_dp)
     292         1222 :       Cx = 1.5_dp*(3.0_dp/4.0_dp/pi)**(1.0_dp/3.0_dp)
     293              : 
     294              : !$OMP     DO
     295              :       DO ii = 1, npoints
     296              :          !! scale densities with 0.5 because code comes from LSD
     297     23239880 :          my_rho = rho(ii)*0.5_dp
     298     23239880 :          my_ndrho = norm_drho(ii)*0.5_dp
     299     23239880 :          IF (my_rho > my_epsilon_rho) THEN
     300     21552092 :             IF (grad_deriv >= 0) THEN
     301     21552092 :                t1 = my_rho**(0.1e1_dp/0.3e1_dp)
     302     21552092 :                t2 = t1*my_rho
     303     21552092 :                t3 = 0.1e1_dp/t2
     304     21552092 :                xx = my_ndrho*MAX(t3, epsilon_rho43)
     305     21552092 :                t5 = my_ndrho**2
     306     21552092 :                t6 = beta*t5
     307     21552092 :                t7 = my_rho**2
     308     21552092 :                t8 = t1**2
     309     21552092 :                t10 = 0.1e1_dp/t8/t7
     310     21552092 :                t11 = beta*my_ndrho
     311     21552092 :                t12 = LOG(xx + SQRT(xx**0.2e1_dp + 0.1e1_dp))
     312     21552092 :                t16 = 0.10e1_dp + 0.60e1_dp*t11*t3*t12
     313     21552092 :                t17 = 0.1e1_dp/t16
     314     21552092 :                t18 = t10*t17
     315     21552092 :                t21 = 0.20e1_dp*Cx + 0.20e1_dp*t6*t18
     316     21552092 :                t22 = SQRT(t21)
     317     21552092 :                t23 = t22*t21
     318     21552092 :                t24 = my_rho*t23
     319     21552092 :                t26 = rootpi
     320     21552092 :                t27 = 0.1e1_dp/t26
     321     21552092 :                t28 = 0.1e1_dp/omega
     322     21552092 :                t29 = 0.1e1_dp/t22
     323     21552092 :                t30 = t28*t29
     324     21552092 :                t31 = t26*t1
     325     21552092 :                t34 = erf(0.3000000000e1_dp*t30*t31)
     326     21552092 :                t36 = omega*t22
     327     21552092 :                t37 = 0.1e1_dp/t1
     328     21552092 :                t38 = t27*t37
     329     21552092 :                t39 = omega**2
     330     21552092 :                t40 = 0.1e1_dp/t39
     331     21552092 :                t41 = 0.1e1_dp/t21
     332     21552092 :                t42 = t40*t41
     333     21552092 :                t43 = pi*t8
     334     21552092 :                t44 = t42*t43
     335     21552092 :                t46 = EXP(-0.8999999998e1_dp*t44)
     336     21552092 :                t47 = t39*t21
     337     21552092 :                t48 = 0.1e1_dp/pi
     338     21552092 :                t49 = 0.1e1_dp/t8
     339     21552092 :                t51 = t46 - 0.10e1_dp
     340     21552092 :                t52 = t48*t49*t51
     341     21552092 :                t55 = t46 - 0.15e1_dp - 0.5555555558e-1_dp*t47*t52
     342     21552092 :                t59 = t26*t34 + 0.3333333334e0_dp*t36*t38*t55
     343     21552092 :                t60 = t27*t59
     344              : 
     345              :                !! Multiply with 2.0 because it code comes from LSD
     346     21552092 :                e_0(ii) = e_0(ii) - 0.2222222224e0_dp*t24*omega*t60*sx*2.0_dp
     347              : 
     348              :             END IF
     349     21552092 :             IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     350     21122817 :                t64 = t23*omega
     351     21122817 :                t68 = my_rho*t22*omega
     352     21122817 :                t69 = t7*my_rho
     353     21122817 :                t71 = 0.1e1_dp/t8/t69
     354     21122817 :                t72 = t71*t17
     355     21122817 :                t75 = t16**2
     356     21122817 :                t76 = 0.1e1_dp/t75
     357     21122817 :                t77 = t10*t76
     358     21122817 :                t79 = 0.1e1_dp/t1/t7
     359     21122817 :                t84 = 1 + t5*t10
     360     21122817 :                t85 = SQRT(t84)
     361     21122817 :                t86 = 0.1e1_dp/t85
     362     21122817 :                t87 = t71*t86
     363     21122817 :                t90 = -0.8000000000e1_dp*t11*t79*t12 - 0.8000000000e1_dp*t6*t87
     364     21122817 :                t91 = t77*t90
     365     21122817 :                t94 = -0.5333333333e1_dp*t6*t72 - 0.20e1_dp*t6*t91
     366     21122817 :                t95 = t60*t94
     367     21122817 :                t98 = omega*t27
     368     21122817 :                t99 = SQRT(0.3141592654e1_dp)
     369     21122817 :                t100 = 0.1e1_dp/t99
     370     21122817 :                t101 = t26*t100
     371     21122817 :                t103 = EXP(-0.9000000000e1_dp*t44)
     372     21122817 :                t104 = 0.1e1_dp/t23
     373     21122817 :                t105 = t28*t104
     374     21122817 :                t109 = t26*t49
     375              :                t112 = -0.1500000000e1_dp*t105*t31*t94 + 0.1000000000e1_dp*t30* &
     376     21122817 :                       t109
     377     21122817 :                t113 = t103*t112
     378     21122817 :                t116 = omega*t29
     379     21122817 :                t117 = t116*t27
     380     21122817 :                t118 = t37*t55
     381     21122817 :                t122 = t27*t3
     382     21122817 :                t126 = t21**2
     383     21122817 :                t127 = 0.1e1_dp/t126
     384     21122817 :                t128 = t40*t127
     385     21122817 :                t130 = t128*t43*t94
     386     21122817 :                t132 = pi*t37
     387     21122817 :                t133 = t42*t132
     388     21122817 :                t135 = 0.8999999998e1_dp*t130 - 0.5999999999e1_dp*t133
     389     21122817 :                t136 = t135*t46
     390     21122817 :                t137 = t39*t94
     391     21122817 :                t140 = t8*my_rho
     392     21122817 :                t141 = 0.1e1_dp/t140
     393     21122817 :                t143 = t48*t141*t51
     394     21122817 :                t146 = t47*t48
     395     21122817 :                t147 = t49*t135
     396     21122817 :                t148 = t147*t46
     397              :                t151 = t136 - 0.5555555558e-1_dp*t137*t52 + 0.3703703705e-1_dp*t47 &
     398     21122817 :                       *t143 - 0.5555555558e-1_dp*t146*t148
     399              :                t155 = REAL(2*t101*t113, KIND=dp) + 0.1666666667e0_dp*t117*t118*t94 - &
     400              :                       0.1111111111e0_dp*t36*t122*t55 + 0.3333333334e0_dp*t36*t38* &
     401     21122817 :                       t151
     402              : 
     403              :                e_rho(ii) = e_rho(ii) + (-0.2222222224e0_dp*t64*t60 - 0.3333333336e0_dp*t68*t95 - &
     404     21122817 :                                         0.2222222224e0_dp*t24*t98*t155)*sx
     405              : 
     406     21122817 :                t168 = 0.60e1_dp*beta*t3*t12 + 0.60e1_dp*t11*t10*t86
     407     21122817 :                t169 = t77*t168
     408     21122817 :                t172 = 0.40e1_dp*t11*t18 - 0.20e1_dp*t6*t169
     409     21122817 :                t173 = t60*t172
     410     21122817 :                t176 = pi*t100
     411     21122817 :                t177 = t176*t103
     412     21122817 :                t185 = t128*pi
     413     21122817 :                t186 = t8*t172
     414     21122817 :                t190 = t39*t172
     415              :                t196 = 0.8999999998e1_dp*t185*t186*t46 - 0.5555555558e-1_dp*t190 &
     416     21122817 :                       *t52 - 0.5000000001e0_dp*t41*t172*t46
     417              :                t200 = -0.3000000000e1_dp*t177*t105*t1*t172 + 0.1666666667e0_dp* &
     418     21122817 :                       t117*t118*t172 + 0.3333333334e0_dp*t36*t38*t196
     419              : 
     420              :                e_ndrho(ii) = e_ndrho(ii) + (-0.3333333336e0_dp*t68*t173 - 0.2222222224e0_dp*t24*t98 &
     421     21122817 :                                             *t200)*sx
     422              : 
     423              :             END IF
     424     21552092 :             IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     425      3660264 :                t207 = t27*t155
     426      3660264 :                t211 = my_rho*t29*omega
     427      3660264 :                t212 = t94**2
     428      3660264 :                t213 = t60*t212
     429      3660264 :                t216 = t207*t94
     430      3660264 :                t219 = t7**2
     431      3660264 :                t221 = 0.1e1_dp/t8/t219
     432      3660264 :                t222 = t221*t17
     433      3660264 :                t225 = t71*t76
     434      3660264 :                t226 = t225*t90
     435      3660264 :                t230 = 0.1e1_dp/t75/t16
     436      3660264 :                t231 = t10*t230
     437      3660264 :                t232 = t90**2
     438      3660264 :                t233 = t231*t232
     439      3660264 :                t237 = 0.1e1_dp/t1/t69
     440      3660264 :                t241 = t221*t86
     441      3660264 :                t244 = t5**2
     442      3660264 :                t245 = beta*t244
     443      3660264 :                t250 = 0.1e1_dp/t85/t84
     444      3660264 :                t251 = 0.1e1_dp/t1/t219/t69*t250
     445              :                t254 = 0.1866666667e2_dp*t11*t237*t12 + 0.4000000000e2_dp*t6*t241 &
     446      3660264 :                       - 0.1066666667e2_dp*t245*t251
     447      3660264 :                t255 = t77*t254
     448              :                t258 = 0.1955555555e2_dp*t6*t222 + 0.1066666667e2_dp*t6*t226 + 0.40e1_dp &
     449      3660264 :                       *t6*t233 - 0.20e1_dp*t6*t255
     450      3660264 :                t259 = t60*t258
     451      3660264 :                t264 = 0.9000000000e1_dp*t130 - 0.6000000000e1_dp*t133
     452      3660264 :                t265 = t264*t103
     453      3660264 :                t270 = 0.1e1_dp/t22/t126
     454      3660264 :                t271 = t28*t270
     455      3660264 :                t281 = t26*t141
     456              :                t284 = 0.2250000000e1_dp*t271*t31*t212 - 0.1000000000e1_dp*t105* &
     457              :                       t109*t94 - 0.1500000000e1_dp*t105*t31*t258 - 0.6666666667e0_dp &
     458      3660264 :                       *t30*t281
     459      3660264 :                t285 = t103*t284
     460      3660264 :                t289 = omega*t104*t27
     461      3660264 :                t293 = t3*t55
     462      3660264 :                t297 = t37*t151
     463      3660264 :                t304 = t27*t79
     464      3660264 :                t311 = t126*t21
     465      3660264 :                t312 = 0.1e1_dp/t311
     466      3660264 :                t313 = t40*t312
     467      3660264 :                t316 = 0.1800000000e2_dp*t313*t43*t212
     468      3660264 :                t319 = 0.1200000000e2_dp*t128*t132*t94
     469      3660264 :                t321 = t128*t43*t258
     470      3660264 :                t323 = pi*t3
     471      3660264 :                t325 = 0.2000000000e1_dp*t42*t323
     472      3660264 :                t326 = -t316 + t319 + 0.8999999998e1_dp*t321 + t325
     473      3660264 :                t328 = t135**2
     474      3660264 :                t330 = t39*t258
     475      3660264 :                t335 = t137*t48
     476      3660264 :                t339 = t48*t10*t51
     477      3660264 :                t343 = t141*t135*t46
     478      3660264 :                t346 = t49*t326
     479      3660264 :                t347 = t346*t46
     480      3660264 :                t351 = t49*t328*t46
     481              :                t354 = t326*t46 + t328*t46 - 0.5555555558e-1_dp*t330*t52 + 0.7407407410e-1_dp &
     482              :                       *t137*t143 - 0.1111111112e0_dp*t335*t148 - 0.6172839508e-1_dp &
     483              :                       *t47*t339 + 0.7407407410e-1_dp*t146*t343 - 0.5555555558e-1_dp &
     484      3660264 :                       *t146*t347 - 0.5555555558e-1_dp*t146*t351
     485              :                t358 = REAL(2*t101*t265*t112, KIND=dp) + REAL(2*t101*t285, KIND=dp) - 0.8333333335e-1_dp &
     486              :                       *t289*t118*t212 - 0.1111111111e0_dp*t117*t293*t94 &
     487              :                       + 0.3333333334e0_dp*t117*t297*t94 + 0.1666666667e0_dp*t117* &
     488              :                       t118*t258 + 0.1481481481e0_dp*t36*t304*t55 - 0.2222222222e0_dp* &
     489      3660264 :                       t36*t122*t151 + 0.3333333334e0_dp*t36*t38*t354
     490              : 
     491              :                e_rho_rho(ii) = e_rho_rho(ii) + (-0.6666666672e0_dp*t36*t95 - 0.4444444448e0_dp*t64*t207 - &
     492              :                                                 0.1666666668e0_dp*t211*t213 - 0.6666666672e0_dp*t68*t216 - 0.3333333336e0_dp &
     493      3660264 :                                                 *t68*t259 - 0.2222222224e0_dp*t24*t98*t358)*sx
     494              : 
     495      3660264 :                t365 = t27*t200
     496      3660264 :                t368 = t94*t172
     497      3660264 :                t372 = t365*t94
     498      3660264 :                t377 = t225*t168
     499      3660264 :                t382 = t6*t10
     500      3660264 :                t383 = t230*t90
     501      3660264 :                t384 = t383*t168
     502      3660264 :                t393 = beta*t5*my_ndrho
     503      3660264 :                t397 = 0.1e1_dp/t1/t219/t7*t250
     504              :                t400 = -0.8000000000e1_dp*beta*t79*t12 - 0.2400000000e2_dp*t11* &
     505      3660264 :                       t87 + 0.8000000000e1_dp*t393*t397
     506      3660264 :                t401 = t77*t400
     507              :                t404 = -0.1066666667e2_dp*t11*t72 + 0.5333333333e1_dp*t6*t377 - 0.40e1_dp &
     508      3660264 :                       *t11*t91 + 0.40e1_dp*t382*t384 - 0.20e1_dp*t6*t401
     509      3660264 :                t405 = t60*t404
     510      3660264 :                t408 = t207*t172
     511      3660264 :                t412 = t26*pi*t100
     512      3660264 :                t413 = t412*t128
     513      3660264 :                t417 = t271*t26
     514      3660264 :                t418 = t1*t94
     515              :                t428 = 0.2250000000e1_dp*t417*t418*t172 - 0.1500000000e1_dp*t105 &
     516      3660264 :                       *t31*t404 - 0.5000000000e0_dp*t105*t109*t172
     517      3660264 :                t429 = t103*t428
     518      3660264 :                t435 = t37*t196
     519      3660264 :                t451 = t313*pi
     520      3660264 :                t452 = t8*t94
     521      3660264 :                t455 = 0.1800000000e2_dp*t451*t452*t172
     522      3660264 :                t457 = t128*t43*t404
     523      3660264 :                t460 = t128*t132*t172
     524      3660264 :                t462 = -t455 + 0.8999999998e1_dp*t457 + 0.5999999999e1_dp*t460
     525      3660264 :                t463 = t462*t46
     526      3660264 :                t464 = t135*t40
     527      3660264 :                t465 = t464*t127
     528      3660264 :                t466 = t172*t46
     529      3660264 :                t467 = t43*t466
     530      3660264 :                t470 = t39*t404
     531      3660264 :                t473 = t94*t127
     532      3660264 :                t478 = 0.1e1_dp/my_rho
     533      3660264 :                t479 = t41*t478
     534      3660264 :                t482 = t190*t48
     535      3660264 :                t486 = t49*t462*t46
     536      3660264 :                t489 = t41*t135
     537              :                t492 = t463 + 0.8999999998e1_dp*t465*t467 - 0.5555555558e-1_dp*t470 &
     538              :                       *t52 - 0.5000000001e0_dp*t473*t466 + 0.3703703705e-1_dp*t190*t143 &
     539              :                       + 0.3333333334e0_dp*t479*t466 - 0.5555555558e-1_dp*t482*t148 &
     540      3660264 :                       - 0.5555555558e-1_dp*t146*t486 - 0.5000000001e0_dp*t489*t466
     541              :                t496 = 0.1800000000e2_dp*t413*t186*t113 + REAL(2*t101*t429, KIND=dp) &
     542              :                       - 0.8333333335e-1_dp*t289*t118*t368 + 0.1666666667e0_dp*t117*t435 &
     543              :                       *t94 + 0.1666666667e0_dp*t117*t118*t404 - 0.5555555555e-1_dp &
     544              :                       *t117*t293*t172 - 0.1111111111e0_dp*t36*t122*t196 + 0.1666666667e0_dp &
     545      3660264 :                       *t117*t297*t172 + 0.3333333334e0_dp*t36*t38*t492
     546              : 
     547              :                e_ndrho_rho(ii) = e_ndrho_rho(ii) + (-0.3333333336e0_dp*t36*t173 - 0.2222222224e0_dp*t64*t365 - &
     548              :                                                     0.1666666668e0_dp*t211*t60*t368 - 0.3333333336e0_dp*t68*t372 &
     549              :                                                     - 0.3333333336e0_dp*t68*t405 - 0.3333333336e0_dp*t68*t408 - 0.2222222224e0_dp &
     550      3660264 :                                                     *t24*t98*t496)*sx
     551              : 
     552      3660264 :                t501 = t172**2
     553      3660264 :                t502 = t60*t501
     554      3660264 :                t505 = t365*t172
     555      3660264 :                t508 = beta*t10
     556      3660264 :                t513 = t168**2
     557      3660264 :                t514 = t231*t513
     558      3660264 :                t519 = t219*my_rho
     559      3660264 :                t521 = 0.1e1_dp/t1/t519
     560      3660264 :                t525 = 0.120e2_dp*t508*t86 - 0.60e1_dp*t6*t521*t250
     561      3660264 :                t526 = t77*t525
     562              :                t529 = 0.40e1_dp*t508*t17 - 0.80e1_dp*t11*t169 + 0.40e1_dp*t6*t514 &
     563      3660264 :                       - 0.20e1_dp*t6*t526
     564      3660264 :                t530 = t60*t529
     565      3660264 :                t533 = pi**2
     566      3660264 :                t534 = t533*t100
     567      3660264 :                t536 = 0.1e1_dp/t39/omega
     568      3660264 :                t537 = t534*t536
     569      3660264 :                t539 = 0.1e1_dp/t22/t311
     570      3660264 :                t562 = t8*t501
     571      3660264 :                t566 = t8*t529
     572      3660264 :                t570 = t39**2
     573      3660264 :                t571 = 0.1e1_dp/t570
     574      3660264 :                t572 = t126**2
     575      3660264 :                t573 = 0.1e1_dp/t572
     576      3660264 :                t574 = t571*t573
     577      3660264 :                t575 = t574*t533
     578      3660264 :                t576 = t2*t501
     579      3660264 :                t580 = t39*t529
     580              :                t586 = -0.2250000000e2_dp*t451*t562*t46 + 0.8999999998e1_dp*t185 &
     581              :                       *t566*t46 + 0.8099999996e2_dp*t575*t576*t46 - 0.5555555558e-1_dp &
     582      3660264 :                       *t580*t52 - 0.5000000001e0_dp*t41*t529*t46
     583              :                t590 = -0.2700000000e2_dp*t537*t539*my_rho*t501*t103 + 0.4500000000e1_dp &
     584              :                       *t177*t271*t1*t501 - 0.3000000000e1_dp*t177*t105* &
     585              :                       t1*t529 - 0.8333333335e-1_dp*t289*t118*t501 + 0.3333333334e0_dp &
     586              :                       *t117*t435*t172 + 0.1666666667e0_dp*t117*t118*t529 + 0.3333333334e0_dp &
     587      3660264 :                       *t36*t38*t586
     588              : 
     589              :                e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + (-0.1666666668e0_dp*t211*t502 - 0.6666666672e0_dp*t68*t505 &
     590              :                                                         - 0.3333333336e0_dp*t68*t530 - 0.2222222224e0_dp*t24*t98*t590) &
     591      3660264 :                                    *sx
     592              : 
     593              :             END IF
     594     21552092 :             IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
     595            0 :                t601 = t27*t358
     596            0 :                t605 = my_rho*t104*omega
     597            0 :                t606 = t212*t94
     598            0 :                t613 = t94*t258
     599            0 :                t624 = 0.1e1_dp/t8/t519
     600            0 :                t628 = t221*t76
     601            0 :                t632 = t71*t230
     602            0 :                t639 = t75**2
     603            0 :                t640 = 0.1e1_dp/t639
     604            0 :                t641 = t10*t640
     605            0 :                t657 = t219**2
     606            0 :                t667 = t84**2
     607            0 :                t669 = 0.1e1_dp/t85/t667
     608              :                t677 = -0.9125925923e2_dp*t6*t624*t17 - 0.5866666667e2_dp*t6*t628 &
     609              :                       *t90 - 0.3200000001e2_dp*t6*t632*t232 + 0.1600000000e2_dp*t6 &
     610              :                       *t225*t254 - 0.120e2_dp*t6*t641*t232*t90 + 0.120e2_dp*t382 &
     611              :                       *t383*t254 - 0.20e1_dp*t6*t77*(-0.6222222223e2_dp*t11/t1/ &
     612              :                                                      t219*t12 - 0.2115555556e3_dp*t6*t624*t86 + 0.1315555556e3_dp* &
     613              :                                                      t245/t1/t657*t250 - 0.4266666668e2_dp*beta*t244*t5/t657 &
     614            0 :                                                      /t69*t669)
     615            0 :                t687 = t28*t539
     616            0 :                t716 = t264**2
     617            0 :                t722 = omega*t270*t27
     618            0 :                t735 = t79*t55
     619            0 :                t739 = t3*t151
     620            0 :                t746 = t37*t354
     621            0 :                t769 = t40*t573
     622              :                t791 = 0.5400000000e2_dp*t769*t43*t606 - 0.3600000000e2_dp*t313* &
     623              :                       t132*t212 - 0.5400000000e2_dp*t451*t452*t258 - 0.6000000000e1_dp &
     624              :                       *t128*t323*t94 + 0.1800000000e2_dp*t128*t132*t258 + 0.8999999998e1_dp &
     625            0 :                       *t128*t43*t677 - 0.2666666667e1_dp*t42*pi*t79
     626            0 :                t793 = t328*t135
     627              :                t838 = REAL(3*t326*t135*t46, KIND=dp) + REAL(t791*t46, KIND=dp) + REAL(t793* &
     628              :                                                           t46, KIND=dp) - 0.5555555558e-1_dp*t39*t677*t52 + 0.1111111112e0_dp*t330 &
     629              :                       *t143 - 0.1666666668e0_dp*t330*t48*t148 - 0.1851851853e0_dp*t137 &
     630              :                       *t339 + 0.2222222223e0_dp*t335*t343 - 0.1666666668e0_dp*t335* &
     631              :                       t347 - 0.1666666668e0_dp*t335*t351 + 0.1646090535e0_dp*t47*t48 &
     632              :                       *t71*t51 - 0.1851851853e0_dp*REAL(t146, KIND=dp)*REAL(t10, KIND=dp)*REAL(t135, KIND=dp) &
     633              :                       *REAL(t46, KIND=dp) + 0.1111111112e0_dp*REAL(t146, KIND=dp)*REAL(t141, KIND=dp)*REAL(t326, KIND=dp) &
     634              :                       *REAL(t46, KIND=dp) + 0.1111111112e0_dp*REAL(t146, KIND=dp)*REAL(t141, KIND=dp)*REAL(t328, KIND=dp) &
     635              :                       *REAL(t46, KIND=dp) - 0.5555555558e-1_dp*REAL(t146, KIND=dp)*REAL(t49, KIND=dp)*REAL(t791, KIND=dp) &
     636              :                       *REAL(t46, KIND=dp) - 0.1666666668e0_dp*REAL(t146, KIND=dp)*REAL(t346, KIND=dp)*REAL(t136, KIND=dp) &
     637              :                       - 0.5555555558e-1_dp*REAL(t146, KIND=dp)*REAL(t49, KIND=dp)*REAL(t793, KIND=dp)* &
     638            0 :                       REAL(t46, KIND=dp)
     639              :                t842 = 0.2e1_dp*t101*(-t316 + t319 + 0.9000000000e1_dp*t321 + t325) &
     640              :                       *t103*t112 + 0.2e1_dp*t101*t103*(-0.5625000000e1_dp*t687*t31 &
     641              :                                                        *t606 + 0.2250000000e1_dp*t271*t109*t212 + 0.6750000000e1_dp* &
     642              :                                                        t417*t418*t258 + 0.1000000000e1_dp*t105*t281*t94 - 0.1500000000e1_dp &
     643              :                                                        *t105*t109*t258 - 0.1500000000e1_dp*t105*t31*t677 &
     644              :                                                        + 0.1111111111e1_dp*t30*t26*t10) + 0.4e1_dp*t101*t265*t284 + &
     645              :                       0.2e1_dp*t101*t716*t103*t112 + 0.1250000000e0_dp*t722*t118 &
     646              :                       *t606 + 0.8333333333e-1_dp*t289*t293*t212 - 0.2500000000e0_dp*t289 &
     647              :                       *t297*t212 - 0.2500000000e0_dp*t289*t118*t613 + 0.2222222222e0_dp &
     648              :                       *t117*t735*t94 - 0.3333333333e0_dp*t117*t739*t94 - &
     649              :                       0.1666666667e0_dp*t117*t293*t258 + 0.5000000001e0_dp*t117*t746 &
     650              :                       *t94 + 0.5000000001e0_dp*t117*t297*t258 + 0.1666666667e0_dp*t117 &
     651              :                       *t118*t677 - 0.3456790122e0_dp*t36*t27*t237*t55 + 0.4444444444e0_dp &
     652              :                       *t36*t304*t151 - 0.3333333333e0_dp*t36*t122*t354 &
     653            0 :                       + 0.3333333334e0_dp*t36*t38*t838
     654              :                t846 = -0.5000000004e0_dp*t116*t213 - 0.2000000001e1_dp*t36*t216 &
     655              :                       - 0.1000000001e1_dp*t36*t259 - 0.6666666672e0_dp*t64*t601 + 0.8333333340e-1_dp &
     656              :                       *t605*t60*t606 - 0.5000000004e0_dp*t211*t207*t212 &
     657              :                       - 0.5000000004e0_dp*t211*t60*t613 - 0.1000000001e1_dp*t68* &
     658              :                       t601*t94 - 0.1000000001e1_dp*t68*t207*t258 - 0.3333333336e0_dp* &
     659            0 :                       t68*t60*t677 - 0.2222222224e0_dp*t24*t98*t842
     660              : 
     661            0 :                e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + t846*sx
     662              : 
     663            0 :                t857 = t27*t496
     664            0 :                t860 = t212*t172
     665            0 :                t867 = t94*t404
     666            0 :                t880 = t258*t172
     667              :                t933 = 0.3911111110e2_dp*t11*t222 - 0.1955555555e2_dp*t6*t628*t168 &
     668              :                       + 0.2133333334e2_dp*t11*t226 - 0.2133333334e2_dp*t6*t71*t384 &
     669              :                       + 0.1066666667e2_dp*t6*t225*t400 + 0.80e1_dp*t11*t233 - 0.120e2_dp &
     670              :                       *t382*t640*t232*t168 + 0.80e1_dp*t382*t383*t400 - 0.40e1_dp &
     671              :                       *t11*t255 + 0.40e1_dp*t382*t230*t254*t168 - 0.20e1_dp* &
     672              :                       t6*t77*(0.1866666667e2_dp*beta*t237*t12 + 0.9866666667e2_dp* &
     673              :                               t11*t241 - 0.8266666668e2_dp*t393*t251 + 0.3200000001e2_dp*beta &
     674            0 :                               *t244*my_ndrho/t657/t7*t669)
     675            0 :                t961 = t687*t26
     676            0 :                t1002 = t3*t196
     677              :                t1009 = 0.2e1_dp*t101*(-t455 + 0.9000000000e1_dp*t457 + 0.6000000000e1_dp &
     678              :                                       *t460)*t103*t112 + 0.1800000000e2_dp*t412*t264*t40*t127 &
     679              :                        *t8*t172*t103*t112 + 0.2e1_dp*t101*t265*t428 + 0.1800000000e2_dp &
     680              :                        *t413*t186*t285 + 0.2e1_dp*t101*t103*(-0.5625000000e1_dp &
     681              :                                                              *t961*t1*t212*t172 + 0.4500000000e1_dp*t417*t418*t404 &
     682              :                                                              + 0.1500000000e1_dp*t417*t49*t94*t172 - 0.1000000000e1_dp* &
     683              :                                                           t105*t109*t404 + 0.2250000000e1_dp*t417*t1*t258*t172 - 0.1500000000e1_dp &
     684              :                                                              *t105*t31*t933 + 0.3333333334e0_dp*t105*t281* &
     685              :                                                              t172) + 0.1250000000e0_dp*t722*t118*t860 - 0.8333333335e-1_dp*t289 &
     686              :                        *t435*t212 - 0.1666666667e0_dp*t289*t118*t867 + 0.5555555555e-1_dp &
     687              :                        *t289*t293*t368 - 0.1111111111e0_dp*t117*t1002*t94 &
     688            0 :                        - 0.1111111111e0_dp*t117*t293*t404
     689            0 :                t1013 = t37*t492
     690            0 :                t1044 = t769*pi
     691              :                t1069 = 0.5400000000e2_dp*t1044*t8*t212*t172 - 0.3600000000e2_dp &
     692              :                        *t451*t452*t404 - 0.2400000000e2_dp*t451*t37*t94*t172 + &
     693              :                        0.1200000000e2_dp*t128*t132*t404 - 0.1800000000e2_dp*t451*t8* &
     694              :                        t258*t172 + 0.8999999998e1_dp*t128*t43*t933 - 0.2000000000e1_dp &
     695            0 :                        *t128*t323*t172
     696            0 :                t1091 = t127*t172*t46
     697              :                t1102 = t1069*t46 + 0.8999999998e1_dp*t326*t40*t127*t467 + REAL(2 &
     698              :                                                                      *t136*t462, KIND=dp) + 0.8999999998e1_dp*t328*t40*t127*t467 - &
     699              :                        0.5555555558e-1_dp*t39*t933*t52 - 0.5000000001e0_dp*t258*t127 &
     700              :                        *t466 + 0.7407407410e-1_dp*t470*t143 + 0.6666666668e0_dp*t94*t478 &
     701              :                        *t1091 - 0.1111111112e0_dp*t470*t48*t148 - 0.1111111112e0_dp &
     702            0 :                        *t335*t486 - 0.1000000001e1_dp*t94*t135*t1091
     703              :                t1136 = -0.6172839508e-1_dp*t190*t339 - 0.5555555556e0_dp*t41/t7 &
     704              :                        *t466 + 0.7407407410e-1_dp*t482*t343 + 0.7407407410e-1_dp*t146* &
     705              :                        t141*t462*t46 + 0.6666666668e0_dp*t479*t135*t172*t46 - 0.5555555558e-1_dp &
     706              :                        *t482*t347 - 0.5555555558e-1_dp*t146*t49*t1069 &
     707              :                        *t46 - 0.5000000001e0_dp*t41*t326*t466 - 0.5555555558e-1_dp*t482 &
     708              :                        *t351 - 0.1111111112e0_dp*t146*t147*t463 - 0.5000000001e0_dp* &
     709            0 :                        t41*t328*t466
     710              :                t1141 = -0.1666666667e0_dp*t289*t297*t368 + 0.3333333334e0_dp*t117 &
     711              :                        *t1013*t94 + 0.3333333334e0_dp*t117*t297*t404 - 0.8333333335e-1_dp &
     712              :                        *t289*t118*t880 + 0.1666666667e0_dp*t117*t435*t258 + &
     713              :                        0.1666666667e0_dp*t117*t118*t933 + 0.7407407405e-1_dp*t117*t735 &
     714              :                        *t172 + 0.1481481481e0_dp*t36*t304*t196 - 0.1111111111e0_dp* &
     715              :                        t117*t739*t172 - 0.2222222222e0_dp*t36*t122*t492 + 0.1666666667e0_dp &
     716              :                        *t117*t746*t172 + 0.3333333334e0_dp*t36*t38*(t1102 &
     717            0 :                                                                     + t1136)
     718              :                t1146 = -0.3333333336e0_dp*t117*t59*t94*t172 - 0.6666666672e0_dp &
     719              :                        *t36*t372 - 0.6666666672e0_dp*t36*t405 - 0.6666666672e0_dp*t36 &
     720              :                        *t408 - 0.4444444448e0_dp*t64*t857 + 0.8333333340e-1_dp*t605*t60 &
     721              :                        *t860 - 0.1666666668e0_dp*t211*t365*t212 - 0.3333333336e0_dp* &
     722              :                        t211*t60*t867 - 0.3333333336e0_dp*t211*t207*t368 - 0.6666666672e0_dp &
     723              :                        *t68*t857*t94 - 0.6666666672e0_dp*t68*t207*t404 - 0.1666666668e0_dp &
     724              :                        *t211*t60*t880 - 0.3333333336e0_dp*t68*t365* &
     725              :                        t258 - 0.3333333336e0_dp*t68*t60*t933 - 0.3333333336e0_dp*t68* &
     726            0 :                        t601*t172 - 0.2222222224e0_dp*t24*t98*(t1009 + t1141)
     727              : 
     728            0 :                e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) + t1146*sx
     729              : 
     730            0 :                t1153 = t27*t590
     731            0 :                t1156 = t94*t501
     732            0 :                t1163 = t404*t172
     733            0 :                t1167 = t94*t529
     734            0 :                t1177 = beta*t71
     735              :                t1220 = -0.1066666667e2_dp*t1177*t17 + 0.2133333334e2_dp*t11*t377 &
     736              :                        - 0.1066666667e2_dp*t6*t632*t513 + 0.5333333333e1_dp*t6*t225 &
     737              :                        *t525 - 0.40e1_dp*t508*t76*t90 + 0.160e2_dp*t11*t10*t384 - &
     738              :                        0.80e1_dp*t11*t401 - 0.120e2_dp*t382*t640*t90*t513 + 0.80e1_dp &
     739              :                        *t382*t230*t400*t168 + 0.40e1_dp*t382*t383*t525 - 0.20e1_dp &
     740              :                        *t6*t77*(-0.3200000000e2_dp*t1177*t86 + 0.4800000000e2_dp*t6 &
     741            0 :                                 *t397 - 0.2400000000e2_dp*t245/t657/my_rho*t669)
     742            0 :                t1284 = t37*t586
     743              :                t1334 = 0.5400000000e2_dp*t1044*t452*t501 - 0.3600000000e2_dp*t451 &
     744              :                        *t8*t404*t172 - 0.1800000000e2_dp*t451*t452*t529 + 0.8999999998e1_dp &
     745              :                        *t128*t43*t1220 - 0.1200000000e2_dp*t313*t132*t501 &
     746            0 :                        + 0.5999999999e1_dp*t128*t132*t529
     747            0 :                t1341 = t501*t46
     748            0 :                t1345 = t529*t46
     749            0 :                t1370 = t40*pi
     750              :                t1396 = t1334*t46 + 0.1800000000e2_dp*t462*t40*t127*t467 - 0.2250000000e2_dp &
     751              :                        *t464*t312*t43*t1341 + 0.8999999998e1_dp*t465 &
     752              :                        *t43*t1345 - 0.1000000000e1_dp*t404*t127*t466 + 0.8099999996e2_dp &
     753              :                        *t135*t571*t573*t533*t2*t1341 - 0.5555555558e-1_dp*t39 &
     754              :                        *t1220*t52 - 0.5000000001e0_dp*t473*t1345 + 0.3333333334e0_dp* &
     755              :                        t479*t1345 + 0.1000000000e1_dp*t94*t312*t1341 - 0.4500000000e1_dp &
     756              :                        *t94*t573*t501*t1370*t8*t46 + 0.3703703705e-1_dp*t580 &
     757              :                        *t143 + 0.3000000000e1_dp*t312*t37*t501*t1370*t46 - 0.5555555558e-1_dp &
     758              :                        *t580*t48*t148 - 0.1111111112e0_dp*t482*t486 - 0.5555555558e-1_dp &
     759              :                        *t146*t49*t1334*t46 - 0.1000000000e1_dp*t41* &
     760            0 :                        t462*t466 - 0.5000000001e0_dp*t489*t1345
     761              :                t1400 = -0.3600000000e2_dp*t412*t313*t562*t113 + 0.1800000000e2_dp &
     762              :                        *t413*t566*t113 + 0.3600000000e2_dp*t413*t186*t429 + 0.1620000000e3_dp &
     763              :                        *t26*t533*t100*t574*t576*t113 + 0.2e1_dp*t101 &
     764              :                        *t103*(-0.5625000000e1_dp*t961*t418*t501 + 0.4500000000e1_dp &
     765              :                               *t417*t1*t404*t172 + 0.2250000000e1_dp*t417*t418*t529 - &
     766              :                               0.1500000000e1_dp*t105*t31*t1220 + 0.7500000000e0_dp*t271*t109 &
     767              :                               *t501 - 0.5000000000e0_dp*t105*t109*t529) + 0.1250000000e0_dp* &
     768              :                        t722*t118*t1156 - 0.1666666667e0_dp*t289*t435*t368 - 0.1666666667e0_dp &
     769              :                        *t289*t118*t1163 - 0.8333333335e-1_dp*t289*t118*t1167 &
     770              :                        + 0.1666666667e0_dp*t117*t1284*t94 + 0.3333333334e0_dp*t117 &
     771              :                        *t435*t404 + 0.1666666667e0_dp*t117*t118*t1220 + 0.2777777778e-1_dp &
     772              :                        *t289*t293*t501 - 0.1111111111e0_dp*t117*t1002*t172 &
     773              :                        - 0.5555555555e-1_dp*t117*t293*t529 - 0.1111111111e0_dp*t36*t122 &
     774              :                        *t586 - 0.8333333335e-1_dp*t289*t297*t501 + 0.3333333334e0_dp &
     775              :                        *t117*t1013*t172 + 0.1666666667e0_dp*t117*t297*t529 + 0.3333333334e0_dp &
     776            0 :                        *t36*t38*t1396
     777              :                t1404 = -0.1666666668e0_dp*t116*t502 - 0.6666666672e0_dp*t36*t505 &
     778              :                        - 0.3333333336e0_dp*t36*t530 - 0.2222222224e0_dp*t64*t1153 + 0.8333333340e-1_dp &
     779              :                        *t605*t60*t1156 - 0.3333333336e0_dp*t211*t365 &
     780              :                        *t368 - 0.3333333336e0_dp*t211*t60*t1163 - 0.1666666668e0_dp*t211 &
     781              :                        *t60*t1167 - 0.3333333336e0_dp*t68*t1153*t94 - 0.6666666672e0_dp &
     782              :                        *t68*t365*t404 - 0.3333333336e0_dp*t68*t60*t1220 - 0.1666666668e0_dp &
     783              :                        *t211*t207*t501 - 0.6666666672e0_dp*t68*t857* &
     784              :                        t172 - 0.3333333336e0_dp*t68*t207*t529 - 0.2222222224e0_dp*t24* &
     785            0 :                        t98*t1400
     786              : 
     787            0 :                e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) + t1404*sx
     788              : 
     789            0 :                t1405 = t501*t172
     790            0 :                t1412 = t172*t529
     791              :                t1449 = -0.120e2_dp*t508*t76*t168 + 0.240e2_dp*t11*t514 - 0.120e2_dp &
     792              :                        *t11*t526 - 0.120e2_dp*t6*t641*t513*t168 + 0.120e2_dp*t382 &
     793              :                        *t230*t168*t525 - 0.20e1_dp*t6*t77*(-0.240e2_dp*beta*t521 &
     794            0 :                                                            *t250*my_ndrho + 0.180e2_dp*t393/t657*t669)
     795            0 :                t1456 = t1405*t103
     796            0 :                t1467 = t533*pi
     797            0 :                t1472 = t572*t21
     798              :                t1553 = 0.1350000000e3_dp*t537/t22/t572*my_rho*t1456 - 0.8100000000e2_dp &
     799              :                        *t534*t536*t539*my_rho*t172*t103*t529 - 0.2430000000e3_dp &
     800              :                        *t1467*t100/t570/omega/t22/t1472*t140*t1456 - &
     801              :                        0.1125000000e2_dp*t177*t687*t1*t1405 + 0.1350000000e2_dp*t176 &
     802              :                        *t103*t28*t270*t1*t1412 - 0.3000000000e1_dp*t177*t105* &
     803              :                        t1*t1449 + 0.1250000000e0_dp*t722*t118*t1405 - 0.2500000000e0_dp &
     804              :                        *t289*t435*t501 - 0.2500000000e0_dp*t289*t118*t1412 + 0.5000000001e0_dp &
     805              :                        *t117*t1284*t172 + 0.5000000001e0_dp*t117*t435 &
     806              :                        *t529 + 0.1666666667e0_dp*t117*t118*t1449 + 0.3333333334e0_dp*t36 &
     807              :                        *t38*(0.6750000000e2_dp*t1044*t8*t1405*t46 - 0.6750000000e2_dp &
     808              :                              *t451*t186*t1345 - 0.5264999998e3_dp*t571/t1472*t533 &
     809              :                              *t2*t1405*t46 + 0.8999999998e1_dp*t185*t8*t1449*t46 + 0.2429999999e3_dp &
     810              :                              *t575*t2*t529*t466 + 0.7289999995e3_dp/t570/t39 &
     811              :                              /t572/t126*t1467*t7*t1405*t46 - 0.5555555558e-1_dp*t39 &
     812            0 :                              *t1449*t52 - 0.5000000001e0_dp*t41*t1449*t46)
     813              : 
     814              :                e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) + (0.8333333340e-1_dp*t605*t60*t1405 - &
     815              :                                            0.5000000004e0_dp*t211*t365*t501 - 0.5000000004e0_dp*t211*t60*t1412 - 0.1000000001e1_dp &
     816              :                                                              *t68*t1153*t172 - 0.1000000001e1_dp*t68*t365*t529 - 0.3333333336e0_dp &
     817              :                                                                     *t68*t60*t1449 - 0.2222222224e0_dp*t24*t98*t1553) &
     818            0 :                                          *sx
     819              : 
     820              :             END IF
     821              :          END IF
     822              :       END DO
     823              : !$OMP     END DO
     824              : 
     825         1222 :    END SUBROUTINE xb88_lr_lda_calc
     826              : 
     827              : ! **************************************************************************************************
     828              : !> \brief evaluates the becke 88 longrange exchange functional for lsd
     829              : !> \param rho_set ...
     830              : !> \param deriv_set ...
     831              : !> \param grad_deriv ...
     832              : !> \param xb88_lr_params ...
     833              : !> \par History
     834              : !>      01.2008 created [mguidon]
     835              : !> \author Manuel Guidon
     836              : ! **************************************************************************************************
     837          280 :    SUBROUTINE xb88_lr_lsd_eval(rho_set, deriv_set, grad_deriv, xb88_lr_params)
     838              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     839              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     840              :       INTEGER, INTENT(in)                                :: grad_deriv
     841              :       TYPE(section_vals_type), POINTER                   :: xb88_lr_params
     842              : 
     843              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xb88_lr_lsd_eval'
     844              : 
     845              :       INTEGER                                            :: handle, i, ispin, npoints
     846              :       INTEGER, DIMENSION(2, 3)                           :: bo
     847              :       REAL(kind=dp)                                      :: epsilon_rho, omega, sx
     848              :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     849           56 :          POINTER                                         :: dummy, e_0
     850          504 :       TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
     851         1008 :          e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, e_rho_rho_rho, &
     852          336 :          norm_drho, rho
     853              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     854              : 
     855           56 :       CALL timeset(routineN, handle)
     856              : 
     857           56 :       CALL cite_reference(Becke1988)
     858              : 
     859           56 :       NULLIFY (deriv)
     860          168 :       DO i = 1, 2
     861          168 :          NULLIFY (norm_drho(i)%array, rho(i)%array)
     862              :       END DO
     863              : 
     864           56 :       CALL section_vals_val_get(xb88_lr_params, "scale_x", r_val=sx)
     865           56 :       CALL section_vals_val_get(xb88_lr_params, "omega", r_val=omega)
     866              : 
     867              :       CALL xc_rho_set_get(rho_set, &
     868              :                           rhoa=rho(1)%array, &
     869              :                           rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
     870              :                           norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
     871           56 :                           local_bounds=bo)
     872           56 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     873              : 
     874              :       ! meaningful default for the arrays we don't need: let us make compiler
     875              :       ! and debugger happy...
     876           56 :       dummy => rho(1)%array
     877              : 
     878           56 :       e_0 => dummy
     879          168 :       DO i = 1, 2
     880          112 :          e_rho(i)%array => dummy
     881          112 :          e_ndrho(i)%array => dummy
     882          112 :          e_rho_rho(i)%array => dummy
     883          112 :          e_ndrho_rho(i)%array => dummy
     884          112 :          e_ndrho_ndrho(i)%array => dummy
     885          112 :          e_rho_rho_rho(i)%array => dummy
     886          112 :          e_ndrho_rho_rho(i)%array => dummy
     887          112 :          e_ndrho_ndrho_rho(i)%array => dummy
     888          168 :          e_ndrho_ndrho_ndrho(i)%array => dummy
     889              :       END DO
     890              : 
     891           56 :       IF (grad_deriv >= 0) THEN
     892              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     893           56 :                                          allocate_deriv=.TRUE.)
     894           56 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     895              :       END IF
     896           56 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     897              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
     898           56 :                                          allocate_deriv=.TRUE.)
     899           56 :          CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
     900              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
     901           56 :                                          allocate_deriv=.TRUE.)
     902           56 :          CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
     903              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
     904           56 :                                          allocate_deriv=.TRUE.)
     905           56 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
     906              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
     907           56 :                                          allocate_deriv=.TRUE.)
     908           56 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
     909              :       END IF
     910           56 :       IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     911              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
     912            0 :                                          allocate_deriv=.TRUE.)
     913            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(1)%array)
     914              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
     915            0 :                                          allocate_deriv=.TRUE.)
     916            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(2)%array)
     917              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
     918            0 :                                          allocate_deriv=.TRUE.)
     919            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(1)%array)
     920              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
     921            0 :                                          allocate_deriv=.TRUE.)
     922            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(2)%array)
     923              :          deriv => xc_dset_get_derivative(deriv_set, &
     924            0 :                                          [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
     925            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(1)%array)
     926              :          deriv => xc_dset_get_derivative(deriv_set, &
     927            0 :                                          [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
     928            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(2)%array)
     929              :       END IF
     930           56 :       IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
     931              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
     932            0 :                                          allocate_deriv=.TRUE.)
     933            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(1)%array)
     934              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
     935            0 :                                          allocate_deriv=.TRUE.)
     936            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(2)%array)
     937              :          deriv => xc_dset_get_derivative(deriv_set, &
     938            0 :                                          [deriv_norm_drhoa, deriv_rhoa, deriv_rhoa], allocate_deriv=.TRUE.)
     939            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(1)%array)
     940              :          deriv => xc_dset_get_derivative(deriv_set, &
     941            0 :                                          [deriv_norm_drhob, deriv_rhob, deriv_rhob], allocate_deriv=.TRUE.)
     942            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(2)%array)
     943              :          deriv => xc_dset_get_derivative(deriv_set, &
     944            0 :                                          [deriv_norm_drhoa, deriv_norm_drhoa, deriv_rhoa], allocate_deriv=.TRUE.)
     945            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(1)%array)
     946              :          deriv => xc_dset_get_derivative(deriv_set, &
     947            0 :                                          [deriv_norm_drhob, deriv_norm_drhob, deriv_rhob], allocate_deriv=.TRUE.)
     948            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(2)%array)
     949              :          deriv => xc_dset_get_derivative(deriv_set, &
     950            0 :                                          [deriv_norm_drhoa, deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
     951            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(1)%array)
     952              :          deriv => xc_dset_get_derivative(deriv_set, &
     953            0 :                                          [deriv_norm_drhob, deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
     954            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(2)%array)
     955              :       END IF
     956           56 :       IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
     957            0 :          CPABORT("derivatives bigger than 3 not implemented")
     958              :       END IF
     959              : 
     960          168 :       DO ispin = 1, 2
     961              : 
     962              : !$OMP        PARALLEL DEFAULT(NONE) &
     963              : !$OMP                 SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
     964              : !$OMP                 SHARED(e_rho_rho, e_ndrho_rho, e_ndrho_ndrho) &
     965              : !$OMP                 SHARED(e_rho_rho_rho, e_ndrho_rho_rho) &
     966              : !$OMP                 SHARED(e_ndrho_ndrho_rho, e_ndrho_ndrho_ndrho) &
     967              : !$OMP                 SHARED(grad_deriv, npoints, epsilon_rho, sx, omega) &
     968          168 : !$OMP                 SHARED(ispin)
     969              : 
     970              :          CALL xb88_lr_lsd_calc( &
     971              :             rho_spin=rho(ispin)%array, &
     972              :             norm_drho_spin=norm_drho(ispin)%array, &
     973              :             e_0=e_0, &
     974              :             e_rho_spin=e_rho(ispin)%array, &
     975              :             e_ndrho_spin=e_ndrho(ispin)%array, &
     976              :             e_rho_rho_spin=e_rho_rho(ispin)%array, &
     977              :             e_ndrho_rho_spin=e_ndrho_rho(ispin)%array, &
     978              :             e_ndrho_ndrho_spin=e_ndrho_ndrho(ispin)%array, &
     979              :             e_rho_rho_rho_spin=e_rho_rho_rho(ispin)%array, &
     980              :             e_ndrho_rho_rho_spin=e_ndrho_rho_rho(ispin)%array, &
     981              :             e_ndrho_ndrho_rho_spin=e_ndrho_ndrho_rho(ispin)%array, &
     982              :             e_ndrho_ndrho_ndrho_spin=e_ndrho_ndrho_ndrho(ispin)%array, &
     983              :             grad_deriv=grad_deriv, npoints=npoints, &
     984              :             epsilon_rho=epsilon_rho, &
     985              :             sx=sx, omega=omega)
     986              : 
     987              : !$OMP        END PARALLEL
     988              : 
     989              :       END DO
     990              : 
     991           56 :       CALL timestop(handle)
     992              : 
     993           56 :    END SUBROUTINE xb88_lr_lsd_eval
     994              : 
     995              : ! **************************************************************************************************
     996              : !> \brief low level calculation of the becke 88 exchange functional for lsd
     997              : !> \param rho_spin alpha or beta spin density
     998              : !> \param norm_drho_spin || grad rho_spin ||
     999              : !> \param e_0 adds to it the local value of the functional
    1000              : !> \param e_rho_spin e_*_spin derivative of the functional wrt. to the variables
    1001              : !>        named where the * is. Everything wrt. to the spin of the arguments.
    1002              : !> \param e_ndrho_spin ...
    1003              : !> \param e_rho_rho_spin ...
    1004              : !> \param e_ndrho_rho_spin ...
    1005              : !> \param e_ndrho_ndrho_spin ...
    1006              : !> \param e_rho_rho_rho_spin ...
    1007              : !> \param e_ndrho_rho_rho_spin ...
    1008              : !> \param e_ndrho_ndrho_rho_spin ...
    1009              : !> \param e_ndrho_ndrho_ndrho_spin ...
    1010              : !> \param grad_deriv ...
    1011              : !> \param npoints ...
    1012              : !> \param epsilon_rho ...
    1013              : !> \param sx scaling-parameter for exchange
    1014              : !> \param omega ...
    1015              : !> \par History
    1016              : !>      01.2008 created [mguidon]
    1017              : !> \author Manuel Guidon
    1018              : !> \note
    1019              : !>      - Derivatives higher than one not tested
    1020              : ! **************************************************************************************************
    1021          112 :    SUBROUTINE xb88_lr_lsd_calc(rho_spin, norm_drho_spin, e_0, &
    1022              :                                e_rho_spin, e_ndrho_spin, e_rho_rho_spin, e_ndrho_rho_spin, &
    1023              :                                e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
    1024              :                                e_ndrho_ndrho_rho_spin, &
    1025              :                                e_ndrho_ndrho_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx, &
    1026              :                                omega)
    1027              :       REAL(kind=dp), DIMENSION(*), INTENT(in)            :: rho_spin, norm_drho_spin
    1028              :       REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin, e_rho_rho_spin, &
    1029              :          e_ndrho_rho_spin, e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
    1030              :          e_ndrho_ndrho_rho_spin, e_ndrho_ndrho_ndrho_spin
    1031              :       INTEGER, INTENT(in)                                :: grad_deriv, npoints
    1032              :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, sx, omega
    1033              : 
    1034              :       INTEGER                                            :: ii
    1035              :       REAL(kind=dp) :: Cx, epsilon_rho43, my_epsilon_rho, ndrho, rho, t1, t10, t100, t1002, t1009, &
    1036              :          t101, t1013, t103, t104, t1044, t105, t1069, t109, t1091, t11, t1102, t112, t113, t1136, &
    1037              :          t1141, t1146, t1153, t1156, t116, t1163, t1167, t117, t1177, t118, t12, t122, t1220, &
    1038              :          t126, t127, t128, t1284, t130, t132, t133, t1334, t1341, t1345, t135, t136, t137, t1370, &
    1039              :          t1396, t140, t1400, t1404, t1405, t141, t1412, t143, t1449, t1456, t146, t1467, t147, &
    1040              :          t1472, t148, t151, t155, t1553, t16, t168, t169, t17, t172, t173, t176, t177, t18, t185, &
    1041              :          t186, t190, t196, t2, t200, t207, t21, t211, t212, t213, t216
    1042              :       REAL(kind=dp) :: t219, t22, t221, t222, t225, t226, t23, t230, t231, t232, t233, t237, t24, &
    1043              :          t241, t244, t245, t250, t251, t254, t255, t258, t259, t26, t264, t265, t27, t270, t271, &
    1044              :          t28, t281, t284, t285, t289, t29, t293, t297, t3, t30, t304, t31, t311, t312, t313, t316, &
    1045              :          t319, t321, t323, t325, t326, t328, t330, t335, t339, t34, t343, t346, t347, t351, t354, &
    1046              :          t358, t36, t365, t368, t37, t372, t377, t38, t382, t383, t384, t39, t393, t397, t40, &
    1047              :          t400, t401, t404, t405, t408, t41, t412, t413, t417, t418, t42, t428, t429, t43, t435, &
    1048              :          t44, t451, t452, t455, t457, t46, t460, t462, t463, t464, t465
    1049              :       REAL(kind=dp) :: t466, t467, t47, t470, t473, t478, t479, t48, t482, t486, t489, t49, t492, &
    1050              :          t496, t5, t501, t502, t505, t508, t51, t513, t514, t519, t52, t521, t525, t526, t529, &
    1051              :          t530, t533, t534, t536, t537, t539, t55, t562, t566, t570, t571, t572, t573, t574, t575, &
    1052              :          t576, t580, t586, t59, t590, t6, t60, t601, t605, t606, t613, t624, t628, t632, t639, &
    1053              :          t64, t640, t641, t657, t667, t669, t677, t68, t687, t69, t7, t71, t716, t72, t722, t735, &
    1054              :          t739, t746, t75, t76, t769, t77, t79, t791, t793, t8, t838, t84, t842, t846, t85, t857, &
    1055              :          t86, t860, t867, t87, t880, t90, t91, t933, t94, t95, t961, t98
    1056              :       REAL(kind=dp) :: t99, xx
    1057              : 
    1058          112 :       my_epsilon_rho = 0.5_dp*epsilon_rho
    1059          112 :       epsilon_rho43 = my_epsilon_rho**(4.0_dp/3.0_dp)
    1060          112 :       Cx = 1.5_dp*(3.0_dp/4.0_dp/pi)**(1.0_dp/3.0_dp)
    1061              : 
    1062          112 : !$OMP     DO
    1063              :       DO ii = 1, npoints
    1064      4548992 :          rho = rho_spin(ii)
    1065      4548992 :          ndrho = norm_drho_spin(ii)
    1066      4548992 :          IF (rho > my_epsilon_rho) THEN
    1067      4540992 :             IF (grad_deriv >= 0) THEN
    1068      4540992 :                t1 = rho**(0.1e1_dp/0.3e1_dp)
    1069      4540992 :                t2 = t1*rho
    1070      4540992 :                t3 = 0.1e1_dp/t2
    1071      4540992 :                xx = ndrho*MAX(t3, epsilon_rho43)
    1072      4540992 :                t5 = ndrho**2
    1073      4540992 :                t6 = beta*t5
    1074      4540992 :                t7 = rho**2
    1075      4540992 :                t8 = t1**2
    1076      4540992 :                t10 = 0.1e1_dp/t8/t7
    1077      4540992 :                t11 = beta*ndrho
    1078      4540992 :                t12 = LOG(xx + SQRT(xx**0.2e1_dp + 0.1e1_dp))
    1079      4540992 :                t16 = 0.10e1_dp + 0.60e1_dp*t11*t3*t12
    1080      4540992 :                t17 = 0.1e1_dp/t16
    1081      4540992 :                t18 = t10*t17
    1082      4540992 :                t21 = 0.20e1_dp*Cx + 0.20e1_dp*t6*t18
    1083      4540992 :                t22 = SQRT(t21)
    1084      4540992 :                t23 = t22*t21
    1085      4540992 :                t24 = rho*t23
    1086      4540992 :                t26 = rootpi
    1087      4540992 :                t27 = 0.1e1_dp/t26
    1088      4540992 :                t28 = 0.1e1_dp/omega
    1089      4540992 :                t29 = 0.1e1_dp/t22
    1090      4540992 :                t30 = t28*t29
    1091      4540992 :                t31 = t26*t1
    1092      4540992 :                t34 = erf(0.3000000000e1_dp*t30*t31)
    1093      4540992 :                t36 = omega*t22
    1094      4540992 :                t37 = 0.1e1_dp/t1
    1095      4540992 :                t38 = t27*t37
    1096      4540992 :                t39 = omega**2
    1097      4540992 :                t40 = 0.1e1_dp/t39
    1098      4540992 :                t41 = 0.1e1_dp/t21
    1099      4540992 :                t42 = t40*t41
    1100      4540992 :                t43 = pi*t8
    1101      4540992 :                t44 = t42*t43
    1102      4540992 :                t46 = EXP(-0.8999999998e1_dp*t44)
    1103      4540992 :                t47 = t39*t21
    1104      4540992 :                t48 = 0.1e1_dp/pi
    1105      4540992 :                t49 = 0.1e1_dp/t8
    1106      4540992 :                t51 = t46 - 0.10e1_dp
    1107      4540992 :                t52 = t48*t49*t51
    1108      4540992 :                t55 = t46 - 0.15e1_dp - 0.5555555558e-1_dp*t47*t52
    1109      4540992 :                t59 = t26*t34 + 0.3333333334e0_dp*t36*t38*t55
    1110      4540992 :                t60 = t27*t59
    1111              : 
    1112      4540992 :                e_0(ii) = e_0(ii) - 0.2222222224e0_dp*t24*omega*t60*sx
    1113              : 
    1114              :             END IF
    1115      4540992 :             IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1116      4540992 :                t64 = t23*omega
    1117      4540992 :                t68 = rho*t22*omega
    1118      4540992 :                t69 = t7*rho
    1119      4540992 :                t71 = 0.1e1_dp/t8/t69
    1120      4540992 :                t72 = t71*t17
    1121      4540992 :                t75 = t16**2
    1122      4540992 :                t76 = 0.1e1_dp/t75
    1123      4540992 :                t77 = t10*t76
    1124      4540992 :                t79 = 0.1e1_dp/t1/t7
    1125      4540992 :                t84 = 1 + t5*t10
    1126      4540992 :                t85 = SQRT(t84)
    1127      4540992 :                t86 = 0.1e1_dp/t85
    1128      4540992 :                t87 = t71*t86
    1129      4540992 :                t90 = -0.8000000000e1_dp*t11*t79*t12 - 0.8000000000e1_dp*t6*t87
    1130      4540992 :                t91 = t77*t90
    1131      4540992 :                t94 = -0.5333333333e1_dp*t6*t72 - 0.20e1_dp*t6*t91
    1132      4540992 :                t95 = t60*t94
    1133      4540992 :                t98 = omega*t27
    1134      4540992 :                t99 = SQRT(0.3141592654e1_dp)
    1135      4540992 :                t100 = 0.1e1_dp/t99
    1136      4540992 :                t101 = t26*t100
    1137      4540992 :                t103 = EXP(-0.9000000000e1_dp*t44)
    1138      4540992 :                t104 = 0.1e1_dp/t23
    1139      4540992 :                t105 = t28*t104
    1140      4540992 :                t109 = t26*t49
    1141              :                t112 = -0.1500000000e1_dp*t105*t31*t94 + 0.1000000000e1_dp*t30* &
    1142      4540992 :                       t109
    1143      4540992 :                t113 = t103*t112
    1144      4540992 :                t116 = omega*t29
    1145      4540992 :                t117 = t116*t27
    1146      4540992 :                t118 = t37*t55
    1147      4540992 :                t122 = t27*t3
    1148      4540992 :                t126 = t21**2
    1149      4540992 :                t127 = 0.1e1_dp/t126
    1150      4540992 :                t128 = t40*t127
    1151      4540992 :                t130 = t128*t43*t94
    1152      4540992 :                t132 = pi*t37
    1153      4540992 :                t133 = t42*t132
    1154      4540992 :                t135 = 0.8999999998e1_dp*t130 - 0.5999999999e1_dp*t133
    1155      4540992 :                t136 = t135*t46
    1156      4540992 :                t137 = t39*t94
    1157      4540992 :                t140 = t8*rho
    1158      4540992 :                t141 = 0.1e1_dp/t140
    1159      4540992 :                t143 = t48*t141*t51
    1160      4540992 :                t146 = t47*t48
    1161      4540992 :                t147 = t49*t135
    1162      4540992 :                t148 = t147*t46
    1163              :                t151 = t136 - 0.5555555558e-1_dp*t137*t52 + 0.3703703705e-1_dp*t47 &
    1164      4540992 :                       *t143 - 0.5555555558e-1_dp*t146*t148
    1165              :                t155 = REAL(2*t101*t113, KIND=dp) + 0.1666666667e0_dp*t117*t118*t94 - &
    1166              :                       0.1111111111e0_dp*t36*t122*t55 + 0.3333333334e0_dp*t36*t38* &
    1167      4540992 :                       t151
    1168              : 
    1169              :                e_rho_spin(ii) = e_rho_spin(ii) + (-0.2222222224e0_dp*t64*t60 - 0.3333333336e0_dp*t68*t95 - &
    1170      4540992 :                                                   0.2222222224e0_dp*t24*t98*t155)*sx
    1171              : 
    1172      4540992 :                t168 = 0.60e1_dp*beta*t3*t12 + 0.60e1_dp*t11*t10*t86
    1173      4540992 :                t169 = t77*t168
    1174      4540992 :                t172 = 0.40e1_dp*t11*t18 - 0.20e1_dp*t6*t169
    1175      4540992 :                t173 = t60*t172
    1176      4540992 :                t176 = pi*t100
    1177      4540992 :                t177 = t176*t103
    1178      4540992 :                t185 = t128*pi
    1179      4540992 :                t186 = t8*t172
    1180      4540992 :                t190 = t39*t172
    1181              :                t196 = 0.8999999998e1_dp*t185*t186*t46 - 0.5555555558e-1_dp*t190 &
    1182      4540992 :                       *t52 - 0.5000000001e0_dp*t41*t172*t46
    1183              :                t200 = -0.3000000000e1_dp*t177*t105*t1*t172 + 0.1666666667e0_dp* &
    1184      4540992 :                       t117*t118*t172 + 0.3333333334e0_dp*t36*t38*t196
    1185              : 
    1186              :                e_ndrho_spin(ii) = e_ndrho_spin(ii) + (-0.3333333336e0_dp*t68*t173 - 0.2222222224e0_dp*t24*t98 &
    1187      4540992 :                                                       *t200)*sx
    1188              : 
    1189              :             END IF
    1190      4540992 :             IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
    1191            0 :                t207 = t27*t155
    1192            0 :                t211 = rho*t29*omega
    1193            0 :                t212 = t94**2
    1194            0 :                t213 = t60*t212
    1195            0 :                t216 = t207*t94
    1196            0 :                t219 = t7**2
    1197            0 :                t221 = 0.1e1_dp/t8/t219
    1198            0 :                t222 = t221*t17
    1199            0 :                t225 = t71*t76
    1200            0 :                t226 = t225*t90
    1201            0 :                t230 = 0.1e1_dp/t75/t16
    1202            0 :                t231 = t10*t230
    1203            0 :                t232 = t90**2
    1204            0 :                t233 = t231*t232
    1205            0 :                t237 = 0.1e1_dp/t1/t69
    1206            0 :                t241 = t221*t86
    1207            0 :                t244 = t5**2
    1208            0 :                t245 = beta*t244
    1209            0 :                t250 = 0.1e1_dp/t85/t84
    1210            0 :                t251 = 0.1e1_dp/t1/t219/t69*t250
    1211              :                t254 = 0.1866666667e2_dp*t11*t237*t12 + 0.4000000000e2_dp*t6*t241 &
    1212            0 :                       - 0.1066666667e2_dp*t245*t251
    1213            0 :                t255 = t77*t254
    1214              :                t258 = 0.1955555555e2_dp*t6*t222 + 0.1066666667e2_dp*t6*t226 + 0.40e1_dp &
    1215            0 :                       *t6*t233 - 0.20e1_dp*t6*t255
    1216            0 :                t259 = t60*t258
    1217            0 :                t264 = 0.9000000000e1_dp*t130 - 0.6000000000e1_dp*t133
    1218            0 :                t265 = t264*t103
    1219            0 :                t270 = 0.1e1_dp/t22/t126
    1220            0 :                t271 = t28*t270
    1221            0 :                t281 = t26*t141
    1222              :                t284 = 0.2250000000e1_dp*t271*t31*t212 - 0.1000000000e1_dp*t105* &
    1223              :                       t109*t94 - 0.1500000000e1_dp*t105*t31*t258 - 0.6666666667e0_dp &
    1224            0 :                       *t30*t281
    1225            0 :                t285 = t103*t284
    1226            0 :                t289 = omega*t104*t27
    1227            0 :                t293 = t3*t55
    1228            0 :                t297 = t37*t151
    1229            0 :                t304 = t27*t79
    1230            0 :                t311 = t126*t21
    1231            0 :                t312 = 0.1e1_dp/t311
    1232            0 :                t313 = t40*t312
    1233            0 :                t316 = 0.1800000000e2_dp*t313*t43*t212
    1234            0 :                t319 = 0.1200000000e2_dp*t128*t132*t94
    1235            0 :                t321 = t128*t43*t258
    1236            0 :                t323 = pi*t3
    1237            0 :                t325 = 0.2000000000e1_dp*t42*t323
    1238            0 :                t326 = -t316 + t319 + 0.8999999998e1_dp*t321 + t325
    1239            0 :                t328 = t135**2
    1240            0 :                t330 = t39*t258
    1241            0 :                t335 = t137*t48
    1242            0 :                t339 = t48*t10*t51
    1243            0 :                t343 = t141*t135*t46
    1244            0 :                t346 = t49*t326
    1245            0 :                t347 = t346*t46
    1246            0 :                t351 = t49*t328*t46
    1247              :                t354 = t326*t46 + t328*t46 - 0.5555555558e-1_dp*t330*t52 + 0.7407407410e-1_dp &
    1248              :                       *t137*t143 - 0.1111111112e0_dp*t335*t148 - 0.6172839508e-1_dp &
    1249              :                       *t47*t339 + 0.7407407410e-1_dp*t146*t343 - 0.5555555558e-1_dp &
    1250            0 :                       *t146*t347 - 0.5555555558e-1_dp*t146*t351
    1251              :                t358 = REAL(2*t101*t265*t112, KIND=dp) + REAL(2*t101*t285, KIND=dp) - 0.8333333335e-1_dp &
    1252              :                       *t289*t118*t212 - 0.1111111111e0_dp*t117*t293*t94 &
    1253              :                       + 0.3333333334e0_dp*t117*t297*t94 + 0.1666666667e0_dp*t117* &
    1254              :                       t118*t258 + 0.1481481481e0_dp*t36*t304*t55 - 0.2222222222e0_dp* &
    1255            0 :                       t36*t122*t151 + 0.3333333334e0_dp*t36*t38*t354
    1256              : 
    1257              :                e_rho_rho_spin(ii) = e_rho_rho_spin(ii) + (-0.6666666672e0_dp*t36*t95 - 0.4444444448e0_dp*t64*t207 - &
    1258              :                                                       0.1666666668e0_dp*t211*t213 - 0.6666666672e0_dp*t68*t216 - 0.3333333336e0_dp &
    1259            0 :                                                           *t68*t259 - 0.2222222224e0_dp*t24*t98*t358)*sx
    1260              : 
    1261            0 :                t365 = t27*t200
    1262            0 :                t368 = t94*t172
    1263            0 :                t372 = t365*t94
    1264            0 :                t377 = t225*t168
    1265            0 :                t382 = t6*t10
    1266            0 :                t383 = t230*t90
    1267            0 :                t384 = t383*t168
    1268            0 :                t393 = beta*t5*ndrho
    1269            0 :                t397 = 0.1e1_dp/t1/t219/t7*t250
    1270              :                t400 = -0.8000000000e1_dp*beta*t79*t12 - 0.2400000000e2_dp*t11* &
    1271            0 :                       t87 + 0.8000000000e1_dp*t393*t397
    1272            0 :                t401 = t77*t400
    1273              :                t404 = -0.1066666667e2_dp*t11*t72 + 0.5333333333e1_dp*t6*t377 - 0.40e1_dp &
    1274            0 :                       *t11*t91 + 0.40e1_dp*t382*t384 - 0.20e1_dp*t6*t401
    1275            0 :                t405 = t60*t404
    1276            0 :                t408 = t207*t172
    1277            0 :                t412 = t26*pi*t100
    1278            0 :                t413 = t412*t128
    1279            0 :                t417 = t271*t26
    1280            0 :                t418 = t1*t94
    1281              :                t428 = 0.2250000000e1_dp*t417*t418*t172 - 0.1500000000e1_dp*t105 &
    1282            0 :                       *t31*t404 - 0.5000000000e0_dp*t105*t109*t172
    1283            0 :                t429 = t103*t428
    1284            0 :                t435 = t37*t196
    1285            0 :                t451 = t313*pi
    1286            0 :                t452 = t8*t94
    1287            0 :                t455 = 0.1800000000e2_dp*t451*t452*t172
    1288            0 :                t457 = t128*t43*t404
    1289            0 :                t460 = t128*t132*t172
    1290            0 :                t462 = -t455 + 0.8999999998e1_dp*t457 + 0.5999999999e1_dp*t460
    1291            0 :                t463 = t462*t46
    1292            0 :                t464 = t135*t40
    1293            0 :                t465 = t464*t127
    1294            0 :                t466 = t172*t46
    1295            0 :                t467 = t43*t466
    1296            0 :                t470 = t39*t404
    1297            0 :                t473 = t94*t127
    1298            0 :                t478 = 0.1e1_dp/rho
    1299            0 :                t479 = t41*t478
    1300            0 :                t482 = t190*t48
    1301            0 :                t486 = t49*t462*t46
    1302            0 :                t489 = t41*t135
    1303              :                t492 = t463 + 0.8999999998e1_dp*t465*t467 - 0.5555555558e-1_dp*t470 &
    1304              :                       *t52 - 0.5000000001e0_dp*t473*t466 + 0.3703703705e-1_dp*t190*t143 &
    1305              :                       + 0.3333333334e0_dp*t479*t466 - 0.5555555558e-1_dp*t482*t148 &
    1306            0 :                       - 0.5555555558e-1_dp*t146*t486 - 0.5000000001e0_dp*t489*t466
    1307              :                t496 = 0.1800000000e2_dp*t413*t186*t113 + REAL(2*t101*t429, KIND=dp) &
    1308              :                       - 0.8333333335e-1_dp*t289*t118*t368 + 0.1666666667e0_dp*t117*t435 &
    1309              :                       *t94 + 0.1666666667e0_dp*t117*t118*t404 - 0.5555555555e-1_dp &
    1310              :                       *t117*t293*t172 - 0.1111111111e0_dp*t36*t122*t196 + 0.1666666667e0_dp &
    1311            0 :                       *t117*t297*t172 + 0.3333333334e0_dp*t36*t38*t492
    1312              : 
    1313              :                e_ndrho_rho_spin(ii) = e_ndrho_rho_spin(ii) + (-0.3333333336e0_dp*t36*t173 - 0.2222222224e0_dp*t64*t365 - &
    1314              :                                                               0.1666666668e0_dp*t211*t60*t368 - 0.3333333336e0_dp*t68*t372 &
    1315              :                                                      - 0.3333333336e0_dp*t68*t405 - 0.3333333336e0_dp*t68*t408 - 0.2222222224e0_dp &
    1316            0 :                                                               *t24*t98*t496)*sx
    1317              : 
    1318            0 :                t501 = t172**2
    1319            0 :                t502 = t60*t501
    1320            0 :                t505 = t365*t172
    1321            0 :                t508 = beta*t10
    1322            0 :                t513 = t168**2
    1323            0 :                t514 = t231*t513
    1324            0 :                t519 = t219*rho
    1325            0 :                t521 = 0.1e1_dp/t1/t519
    1326            0 :                t525 = 0.120e2_dp*t508*t86 - 0.60e1_dp*t6*t521*t250
    1327            0 :                t526 = t77*t525
    1328              :                t529 = 0.40e1_dp*t508*t17 - 0.80e1_dp*t11*t169 + 0.40e1_dp*t6*t514 &
    1329            0 :                       - 0.20e1_dp*t6*t526
    1330            0 :                t530 = t60*t529
    1331            0 :                t533 = pi**2
    1332            0 :                t534 = t533*t100
    1333            0 :                t536 = 0.1e1_dp/t39/omega
    1334            0 :                t537 = t534*t536
    1335            0 :                t539 = 0.1e1_dp/t22/t311
    1336            0 :                t562 = t8*t501
    1337            0 :                t566 = t8*t529
    1338            0 :                t570 = t39**2
    1339            0 :                t571 = 0.1e1_dp/t570
    1340            0 :                t572 = t126**2
    1341            0 :                t573 = 0.1e1_dp/t572
    1342            0 :                t574 = t571*t573
    1343            0 :                t575 = t574*t533
    1344            0 :                t576 = t2*t501
    1345            0 :                t580 = t39*t529
    1346              :                t586 = -0.2250000000e2_dp*t451*t562*t46 + 0.8999999998e1_dp*t185 &
    1347              :                       *t566*t46 + 0.8099999996e2_dp*t575*t576*t46 - 0.5555555558e-1_dp &
    1348            0 :                       *t580*t52 - 0.5000000001e0_dp*t41*t529*t46
    1349              :                t590 = -0.2700000000e2_dp*t537*t539*rho*t501*t103 + 0.4500000000e1_dp &
    1350              :                       *t177*t271*t1*t501 - 0.3000000000e1_dp*t177*t105* &
    1351              :                       t1*t529 - 0.8333333335e-1_dp*t289*t118*t501 + 0.3333333334e0_dp &
    1352              :                       *t117*t435*t172 + 0.1666666667e0_dp*t117*t118*t529 + 0.3333333334e0_dp &
    1353            0 :                       *t36*t38*t586
    1354              : 
    1355              :                e_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_spin(ii) + (-0.1666666668e0_dp*t211*t502 - 0.6666666672e0_dp*t68*t505 &
    1356              :                                                                   - 0.3333333336e0_dp*t68*t530 - 0.2222222224e0_dp*t24*t98*t590) &
    1357            0 :                                         *sx
    1358              : 
    1359              :             END IF
    1360      4540992 :             IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
    1361            0 :                t601 = t27*t358
    1362            0 :                t605 = rho*t104*omega
    1363            0 :                t606 = t212*t94
    1364            0 :                t613 = t94*t258
    1365            0 :                t624 = 0.1e1_dp/t8/t519
    1366            0 :                t628 = t221*t76
    1367            0 :                t632 = t71*t230
    1368            0 :                t639 = t75**2
    1369            0 :                t640 = 0.1e1_dp/t639
    1370            0 :                t641 = t10*t640
    1371            0 :                t657 = t219**2
    1372            0 :                t667 = t84**2
    1373            0 :                t669 = 0.1e1_dp/t85/t667
    1374              :                t677 = -0.9125925923e2_dp*t6*t624*t17 - 0.5866666667e2_dp*t6*t628 &
    1375              :                       *t90 - 0.3200000001e2_dp*t6*t632*t232 + 0.1600000000e2_dp*t6 &
    1376              :                       *t225*t254 - 0.120e2_dp*t6*t641*t232*t90 + 0.120e2_dp*t382 &
    1377              :                       *t383*t254 - 0.20e1_dp*t6*t77*(-0.6222222223e2_dp*t11/t1/ &
    1378              :                                                      t219*t12 - 0.2115555556e3_dp*t6*t624*t86 + 0.1315555556e3_dp* &
    1379              :                                                      t245/t1/t657*t250 - 0.4266666668e2_dp*beta*t244*t5/t657 &
    1380            0 :                                                      /t69*t669)
    1381            0 :                t687 = t28*t539
    1382            0 :                t716 = t264**2
    1383            0 :                t722 = omega*t270*t27
    1384            0 :                t735 = t79*t55
    1385            0 :                t739 = t3*t151
    1386            0 :                t746 = t37*t354
    1387            0 :                t769 = t40*t573
    1388              :                t791 = 0.5400000000e2_dp*t769*t43*t606 - 0.3600000000e2_dp*t313* &
    1389              :                       t132*t212 - 0.5400000000e2_dp*t451*t452*t258 - 0.6000000000e1_dp &
    1390              :                       *t128*t323*t94 + 0.1800000000e2_dp*t128*t132*t258 + 0.8999999998e1_dp &
    1391            0 :                       *t128*t43*t677 - 0.2666666667e1_dp*t42*pi*t79
    1392            0 :                t793 = t328*t135
    1393              :                t838 = REAL(3*t326*t135*t46, KIND=dp) + REAL(t791*t46, KIND=dp) + REAL(t793* &
    1394              :                                                           t46, KIND=dp) - 0.5555555558e-1_dp*t39*t677*t52 + 0.1111111112e0_dp*t330 &
    1395              :                       *t143 - 0.1666666668e0_dp*t330*t48*t148 - 0.1851851853e0_dp*t137 &
    1396              :                       *t339 + 0.2222222223e0_dp*t335*t343 - 0.1666666668e0_dp*t335* &
    1397              :                       t347 - 0.1666666668e0_dp*t335*t351 + 0.1646090535e0_dp*t47*t48 &
    1398              :                       *t71*t51 - 0.1851851853e0_dp*REAL(t146, KIND=dp)*REAL(t10, KIND=dp)*REAL(t135, KIND=dp) &
    1399              :                       *REAL(t46, KIND=dp) + 0.1111111112e0_dp*REAL(t146, KIND=dp)*REAL(t141, KIND=dp)*REAL(t326, KIND=dp) &
    1400              :                       *REAL(t46, KIND=dp) + 0.1111111112e0_dp*REAL(t146, KIND=dp)*REAL(t141, KIND=dp)*REAL(t328, KIND=dp) &
    1401              :                       *REAL(t46, KIND=dp) - 0.5555555558e-1_dp*REAL(t146, KIND=dp)*REAL(t49, KIND=dp)*REAL(t791, KIND=dp) &
    1402              :                       *REAL(t46, KIND=dp) - 0.1666666668e0_dp*REAL(t146, KIND=dp)*REAL(t346, KIND=dp)*REAL(t136, KIND=dp) &
    1403              :                       - 0.5555555558e-1_dp*REAL(t146, KIND=dp)*REAL(t49, KIND=dp)*REAL(t793, KIND=dp)* &
    1404            0 :                       REAL(t46, KIND=dp)
    1405              :                t842 = 0.2e1_dp*t101*(-t316 + t319 + 0.9000000000e1_dp*t321 + t325) &
    1406              :                       *t103*t112 + 0.2e1_dp*t101*t103*(-0.5625000000e1_dp*t687*t31 &
    1407              :                                                        *t606 + 0.2250000000e1_dp*t271*t109*t212 + 0.6750000000e1_dp* &
    1408              :                                                        t417*t418*t258 + 0.1000000000e1_dp*t105*t281*t94 - 0.1500000000e1_dp &
    1409              :                                                        *t105*t109*t258 - 0.1500000000e1_dp*t105*t31*t677 &
    1410              :                                                        + 0.1111111111e1_dp*t30*t26*t10) + 0.4e1_dp*t101*t265*t284 + &
    1411              :                       0.2e1_dp*t101*t716*t103*t112 + 0.1250000000e0_dp*t722*t118 &
    1412              :                       *t606 + 0.8333333333e-1_dp*t289*t293*t212 - 0.2500000000e0_dp*t289 &
    1413              :                       *t297*t212 - 0.2500000000e0_dp*t289*t118*t613 + 0.2222222222e0_dp &
    1414              :                       *t117*t735*t94 - 0.3333333333e0_dp*t117*t739*t94 - &
    1415              :                       0.1666666667e0_dp*t117*t293*t258 + 0.5000000001e0_dp*t117*t746 &
    1416              :                       *t94 + 0.5000000001e0_dp*t117*t297*t258 + 0.1666666667e0_dp*t117 &
    1417              :                       *t118*t677 - 0.3456790122e0_dp*t36*t27*t237*t55 + 0.4444444444e0_dp &
    1418              :                       *t36*t304*t151 - 0.3333333333e0_dp*t36*t122*t354 &
    1419            0 :                       + 0.3333333334e0_dp*t36*t38*t838
    1420              :                t846 = -0.5000000004e0_dp*t116*t213 - 0.2000000001e1_dp*t36*t216 &
    1421              :                       - 0.1000000001e1_dp*t36*t259 - 0.6666666672e0_dp*t64*t601 + 0.8333333340e-1_dp &
    1422              :                       *t605*t60*t606 - 0.5000000004e0_dp*t211*t207*t212 &
    1423              :                       - 0.5000000004e0_dp*t211*t60*t613 - 0.1000000001e1_dp*t68* &
    1424              :                       t601*t94 - 0.1000000001e1_dp*t68*t207*t258 - 0.3333333336e0_dp* &
    1425            0 :                       t68*t60*t677 - 0.2222222224e0_dp*t24*t98*t842
    1426              : 
    1427            0 :                e_rho_rho_rho_spin(ii) = e_rho_rho_rho_spin(ii) + t846*sx
    1428              : 
    1429            0 :                t857 = t27*t496
    1430            0 :                t860 = t212*t172
    1431            0 :                t867 = t94*t404
    1432            0 :                t880 = t258*t172
    1433              :                t933 = 0.3911111110e2_dp*t11*t222 - 0.1955555555e2_dp*t6*t628*t168 &
    1434              :                       + 0.2133333334e2_dp*t11*t226 - 0.2133333334e2_dp*t6*t71*t384 &
    1435              :                       + 0.1066666667e2_dp*t6*t225*t400 + 0.80e1_dp*t11*t233 - 0.120e2_dp &
    1436              :                       *t382*t640*t232*t168 + 0.80e1_dp*t382*t383*t400 - 0.40e1_dp &
    1437              :                       *t11*t255 + 0.40e1_dp*t382*t230*t254*t168 - 0.20e1_dp* &
    1438              :                       t6*t77*(0.1866666667e2_dp*beta*t237*t12 + 0.9866666667e2_dp* &
    1439              :                               t11*t241 - 0.8266666668e2_dp*t393*t251 + 0.3200000001e2_dp*beta &
    1440            0 :                               *t244*ndrho/t657/t7*t669)
    1441            0 :                t961 = t687*t26
    1442            0 :                t1002 = t3*t196
    1443              :                t1009 = 0.2e1_dp*t101*(-t455 + 0.9000000000e1_dp*t457 + 0.6000000000e1_dp &
    1444              :                                       *t460)*t103*t112 + 0.1800000000e2_dp*t412*t264*t40*t127 &
    1445              :                        *t8*t172*t103*t112 + 0.2e1_dp*t101*t265*t428 + 0.1800000000e2_dp &
    1446              :                        *t413*t186*t285 + 0.2e1_dp*t101*t103*(-0.5625000000e1_dp &
    1447              :                                                              *t961*t1*t212*t172 + 0.4500000000e1_dp*t417*t418*t404 &
    1448              :                                                              + 0.1500000000e1_dp*t417*t49*t94*t172 - 0.1000000000e1_dp* &
    1449              :                                                           t105*t109*t404 + 0.2250000000e1_dp*t417*t1*t258*t172 - 0.1500000000e1_dp &
    1450              :                                                              *t105*t31*t933 + 0.3333333334e0_dp*t105*t281* &
    1451              :                                                              t172) + 0.1250000000e0_dp*t722*t118*t860 - 0.8333333335e-1_dp*t289 &
    1452              :                        *t435*t212 - 0.1666666667e0_dp*t289*t118*t867 + 0.5555555555e-1_dp &
    1453              :                        *t289*t293*t368 - 0.1111111111e0_dp*t117*t1002*t94 &
    1454            0 :                        - 0.1111111111e0_dp*t117*t293*t404
    1455            0 :                t1013 = t37*t492
    1456            0 :                t1044 = t769*pi
    1457              :                t1069 = 0.5400000000e2_dp*t1044*t8*t212*t172 - 0.3600000000e2_dp &
    1458              :                        *t451*t452*t404 - 0.2400000000e2_dp*t451*t37*t94*t172 + &
    1459              :                        0.1200000000e2_dp*t128*t132*t404 - 0.1800000000e2_dp*t451*t8* &
    1460              :                        t258*t172 + 0.8999999998e1_dp*t128*t43*t933 - 0.2000000000e1_dp &
    1461            0 :                        *t128*t323*t172
    1462            0 :                t1091 = t127*t172*t46
    1463              :                t1102 = t1069*t46 + 0.8999999998e1_dp*t326*t40*t127*t467 + REAL(2 &
    1464              :                                                                      *t136*t462, KIND=dp) + 0.8999999998e1_dp*t328*t40*t127*t467 - &
    1465              :                        0.5555555558e-1_dp*t39*t933*t52 - 0.5000000001e0_dp*t258*t127 &
    1466              :                        *t466 + 0.7407407410e-1_dp*t470*t143 + 0.6666666668e0_dp*t94*t478 &
    1467              :                        *t1091 - 0.1111111112e0_dp*t470*t48*t148 - 0.1111111112e0_dp &
    1468            0 :                        *t335*t486 - 0.1000000001e1_dp*t94*t135*t1091
    1469              :                t1136 = -0.6172839508e-1_dp*t190*t339 - 0.5555555556e0_dp*t41/t7 &
    1470              :                        *t466 + 0.7407407410e-1_dp*t482*t343 + 0.7407407410e-1_dp*t146* &
    1471              :                        t141*t462*t46 + 0.6666666668e0_dp*t479*t135*t172*t46 - 0.5555555558e-1_dp &
    1472              :                        *t482*t347 - 0.5555555558e-1_dp*t146*t49*t1069 &
    1473              :                        *t46 - 0.5000000001e0_dp*t41*t326*t466 - 0.5555555558e-1_dp*t482 &
    1474              :                        *t351 - 0.1111111112e0_dp*t146*t147*t463 - 0.5000000001e0_dp* &
    1475            0 :                        t41*t328*t466
    1476              :                t1141 = -0.1666666667e0_dp*t289*t297*t368 + 0.3333333334e0_dp*t117 &
    1477              :                        *t1013*t94 + 0.3333333334e0_dp*t117*t297*t404 - 0.8333333335e-1_dp &
    1478              :                        *t289*t118*t880 + 0.1666666667e0_dp*t117*t435*t258 + &
    1479              :                        0.1666666667e0_dp*t117*t118*t933 + 0.7407407405e-1_dp*t117*t735 &
    1480              :                        *t172 + 0.1481481481e0_dp*t36*t304*t196 - 0.1111111111e0_dp* &
    1481              :                        t117*t739*t172 - 0.2222222222e0_dp*t36*t122*t492 + 0.1666666667e0_dp &
    1482              :                        *t117*t746*t172 + 0.3333333334e0_dp*t36*t38*(t1102 &
    1483            0 :                                                                     + t1136)
    1484              :                t1146 = -0.3333333336e0_dp*t117*t59*t94*t172 - 0.6666666672e0_dp &
    1485              :                        *t36*t372 - 0.6666666672e0_dp*t36*t405 - 0.6666666672e0_dp*t36 &
    1486              :                        *t408 - 0.4444444448e0_dp*t64*t857 + 0.8333333340e-1_dp*t605*t60 &
    1487              :                        *t860 - 0.1666666668e0_dp*t211*t365*t212 - 0.3333333336e0_dp* &
    1488              :                        t211*t60*t867 - 0.3333333336e0_dp*t211*t207*t368 - 0.6666666672e0_dp &
    1489              :                        *t68*t857*t94 - 0.6666666672e0_dp*t68*t207*t404 - 0.1666666668e0_dp &
    1490              :                        *t211*t60*t880 - 0.3333333336e0_dp*t68*t365* &
    1491              :                        t258 - 0.3333333336e0_dp*t68*t60*t933 - 0.3333333336e0_dp*t68* &
    1492            0 :                        t601*t172 - 0.2222222224e0_dp*t24*t98*(t1009 + t1141)
    1493              : 
    1494            0 :                e_ndrho_rho_rho_spin(ii) = e_ndrho_rho_rho_spin(ii) + t1146*sx
    1495              : 
    1496            0 :                t1153 = t27*t590
    1497            0 :                t1156 = t94*t501
    1498            0 :                t1163 = t404*t172
    1499            0 :                t1167 = t94*t529
    1500            0 :                t1177 = beta*t71
    1501              :                t1220 = -0.1066666667e2_dp*t1177*t17 + 0.2133333334e2_dp*t11*t377 &
    1502              :                        - 0.1066666667e2_dp*t6*t632*t513 + 0.5333333333e1_dp*t6*t225 &
    1503              :                        *t525 - 0.40e1_dp*t508*t76*t90 + 0.160e2_dp*t11*t10*t384 - &
    1504              :                        0.80e1_dp*t11*t401 - 0.120e2_dp*t382*t640*t90*t513 + 0.80e1_dp &
    1505              :                        *t382*t230*t400*t168 + 0.40e1_dp*t382*t383*t525 - 0.20e1_dp &
    1506              :                        *t6*t77*(-0.3200000000e2_dp*t1177*t86 + 0.4800000000e2_dp*t6 &
    1507            0 :                                 *t397 - 0.2400000000e2_dp*t245/t657/rho*t669)
    1508            0 :                t1284 = t37*t586
    1509              :                t1334 = 0.5400000000e2_dp*t1044*t452*t501 - 0.3600000000e2_dp*t451 &
    1510              :                        *t8*t404*t172 - 0.1800000000e2_dp*t451*t452*t529 + 0.8999999998e1_dp &
    1511              :                        *t128*t43*t1220 - 0.1200000000e2_dp*t313*t132*t501 &
    1512            0 :                        + 0.5999999999e1_dp*t128*t132*t529
    1513            0 :                t1341 = t501*t46
    1514            0 :                t1345 = t529*t46
    1515            0 :                t1370 = t40*pi
    1516              :                t1396 = t1334*t46 + 0.1800000000e2_dp*t462*t40*t127*t467 - 0.2250000000e2_dp &
    1517              :                        *t464*t312*t43*t1341 + 0.8999999998e1_dp*t465 &
    1518              :                        *t43*t1345 - 0.1000000000e1_dp*t404*t127*t466 + 0.8099999996e2_dp &
    1519              :                        *t135*t571*t573*t533*t2*t1341 - 0.5555555558e-1_dp*t39 &
    1520              :                        *t1220*t52 - 0.5000000001e0_dp*t473*t1345 + 0.3333333334e0_dp* &
    1521              :                        t479*t1345 + 0.1000000000e1_dp*t94*t312*t1341 - 0.4500000000e1_dp &
    1522              :                        *t94*t573*t501*t1370*t8*t46 + 0.3703703705e-1_dp*t580 &
    1523              :                        *t143 + 0.3000000000e1_dp*t312*t37*t501*t1370*t46 - 0.5555555558e-1_dp &
    1524              :                        *t580*t48*t148 - 0.1111111112e0_dp*t482*t486 - 0.5555555558e-1_dp &
    1525              :                        *t146*t49*t1334*t46 - 0.1000000000e1_dp*t41* &
    1526            0 :                        t462*t466 - 0.5000000001e0_dp*t489*t1345
    1527              :                t1400 = -0.3600000000e2_dp*t412*t313*t562*t113 + 0.1800000000e2_dp &
    1528              :                        *t413*t566*t113 + 0.3600000000e2_dp*t413*t186*t429 + 0.1620000000e3_dp &
    1529              :                        *t26*t533*t100*t574*t576*t113 + 0.2e1_dp*t101 &
    1530              :                        *t103*(-0.5625000000e1_dp*t961*t418*t501 + 0.4500000000e1_dp &
    1531              :                               *t417*t1*t404*t172 + 0.2250000000e1_dp*t417*t418*t529 - &
    1532              :                               0.1500000000e1_dp*t105*t31*t1220 + 0.7500000000e0_dp*t271*t109 &
    1533              :                               *t501 - 0.5000000000e0_dp*t105*t109*t529) + 0.1250000000e0_dp* &
    1534              :                        t722*t118*t1156 - 0.1666666667e0_dp*t289*t435*t368 - 0.1666666667e0_dp &
    1535              :                        *t289*t118*t1163 - 0.8333333335e-1_dp*t289*t118*t1167 &
    1536              :                        + 0.1666666667e0_dp*t117*t1284*t94 + 0.3333333334e0_dp*t117 &
    1537              :                        *t435*t404 + 0.1666666667e0_dp*t117*t118*t1220 + 0.2777777778e-1_dp &
    1538              :                        *t289*t293*t501 - 0.1111111111e0_dp*t117*t1002*t172 &
    1539              :                        - 0.5555555555e-1_dp*t117*t293*t529 - 0.1111111111e0_dp*t36*t122 &
    1540              :                        *t586 - 0.8333333335e-1_dp*t289*t297*t501 + 0.3333333334e0_dp &
    1541              :                        *t117*t1013*t172 + 0.1666666667e0_dp*t117*t297*t529 + 0.3333333334e0_dp &
    1542            0 :                        *t36*t38*t1396
    1543              :                t1404 = -0.1666666668e0_dp*t116*t502 - 0.6666666672e0_dp*t36*t505 &
    1544              :                        - 0.3333333336e0_dp*t36*t530 - 0.2222222224e0_dp*t64*t1153 + 0.8333333340e-1_dp &
    1545              :                        *t605*t60*t1156 - 0.3333333336e0_dp*t211*t365 &
    1546              :                        *t368 - 0.3333333336e0_dp*t211*t60*t1163 - 0.1666666668e0_dp*t211 &
    1547              :                        *t60*t1167 - 0.3333333336e0_dp*t68*t1153*t94 - 0.6666666672e0_dp &
    1548              :                        *t68*t365*t404 - 0.3333333336e0_dp*t68*t60*t1220 - 0.1666666668e0_dp &
    1549              :                        *t211*t207*t501 - 0.6666666672e0_dp*t68*t857* &
    1550              :                        t172 - 0.3333333336e0_dp*t68*t207*t529 - 0.2222222224e0_dp*t24* &
    1551            0 :                        t98*t1400
    1552              : 
    1553            0 :                e_ndrho_ndrho_rho_spin(ii) = e_ndrho_ndrho_rho_spin(ii) + t1404*sx
    1554              : 
    1555            0 :                t1405 = t501*t172
    1556            0 :                t1412 = t172*t529
    1557              :                t1449 = -0.120e2_dp*t508*t76*t168 + 0.240e2_dp*t11*t514 - 0.120e2_dp &
    1558              :                        *t11*t526 - 0.120e2_dp*t6*t641*t513*t168 + 0.120e2_dp*t382 &
    1559              :                        *t230*t168*t525 - 0.20e1_dp*t6*t77*(-0.240e2_dp*beta*t521 &
    1560            0 :                                                            *t250*ndrho + 0.180e2_dp*t393/t657*t669)
    1561            0 :                t1456 = t1405*t103
    1562            0 :                t1467 = t533*pi
    1563            0 :                t1472 = t572*t21
    1564              :                t1553 = 0.1350000000e3_dp*t537/t22/t572*rho*t1456 - 0.8100000000e2_dp &
    1565              :                        *t534*t536*t539*rho*t172*t103*t529 - 0.2430000000e3_dp &
    1566              :                        *t1467*t100/t570/omega/t22/t1472*t140*t1456 - &
    1567              :                        0.1125000000e2_dp*t177*t687*t1*t1405 + 0.1350000000e2_dp*t176 &
    1568              :                        *t103*t28*t270*t1*t1412 - 0.3000000000e1_dp*t177*t105* &
    1569              :                        t1*t1449 + 0.1250000000e0_dp*t722*t118*t1405 - 0.2500000000e0_dp &
    1570              :                        *t289*t435*t501 - 0.2500000000e0_dp*t289*t118*t1412 + 0.5000000001e0_dp &
    1571              :                        *t117*t1284*t172 + 0.5000000001e0_dp*t117*t435 &
    1572              :                        *t529 + 0.1666666667e0_dp*t117*t118*t1449 + 0.3333333334e0_dp*t36 &
    1573              :                        *t38*(0.6750000000e2_dp*t1044*t8*t1405*t46 - 0.6750000000e2_dp &
    1574              :                              *t451*t186*t1345 - 0.5264999998e3_dp*t571/t1472*t533 &
    1575              :                              *t2*t1405*t46 + 0.8999999998e1_dp*t185*t8*t1449*t46 + 0.2429999999e3_dp &
    1576              :                              *t575*t2*t529*t466 + 0.7289999995e3_dp/t570/t39 &
    1577              :                              /t572/t126*t1467*t7*t1405*t46 - 0.5555555558e-1_dp*t39 &
    1578            0 :                              *t1449*t52 - 0.5000000001e0_dp*t41*t1449*t46)
    1579              : 
    1580              :                e_ndrho_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_ndrho_spin(ii) + (0.8333333340e-1_dp*t605*t60*t1405 - &
    1581              :                                            0.5000000004e0_dp*t211*t365*t501 - 0.5000000004e0_dp*t211*t60*t1412 - 0.1000000001e1_dp &
    1582              :                                                              *t68*t1153*t172 - 0.1000000001e1_dp*t68*t365*t529 - 0.3333333336e0_dp &
    1583              :                                                                               *t68*t60*t1449 - 0.2222222224e0_dp*t24*t98*t1553) &
    1584            0 :                                               *sx
    1585              : 
    1586              :             END IF
    1587              :          END IF
    1588              :       END DO
    1589              : !$OMP     END DO
    1590              : 
    1591          112 :    END SUBROUTINE xb88_lr_lsd_calc
    1592              : 
    1593              : END MODULE xc_xbecke88_long_range
        

Generated by: LCOV version 2.0-1