LCOV - code coverage report
Current view: top level - src/xc - xc_xbecke88_long_range.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 472 771 61.2 %
Date: 2024-03-28 07:31:50 Functions: 6 6 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief calculates the 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        1298 :    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        1298 :       IF (PRESENT(reference)) THEN
      65           8 :          reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LDA version}"
      66             :       END IF
      67        1298 :       IF (PRESENT(shortform)) THEN
      68           8 :          shortform = "Becke 1988 Exchange Functional (LDA)"
      69             :       END IF
      70        1298 :       IF (PRESENT(needs)) THEN
      71        1290 :          needs%rho = .TRUE.
      72        1290 :          needs%norm_drho = .TRUE.
      73             :       END IF
      74        1298 :       IF (PRESENT(max_deriv)) max_deriv = 3
      75             : 
      76        1298 :    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        4896 :    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        1224 :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
     134        1224 :          e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
     135        1224 :          e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
     136             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     137             : 
     138        1224 :       CALL timeset(routineN, handle)
     139             : 
     140        1224 :       CALL section_vals_val_get(xb88_lr_params, "scale_x", r_val=sx)
     141        1224 :       CALL section_vals_val_get(xb88_lr_params, "omega", r_val=omega)
     142             : 
     143        1224 :       CALL cite_reference(Becke1988)
     144             : 
     145             :       CALL xc_rho_set_get(rho_set, rho=rho, &
     146        1224 :                           norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
     147        1224 :       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        1224 :       dummy => rho
     152             : 
     153        1224 :       e_0 => dummy
     154        1224 :       e_rho => dummy
     155        1224 :       e_ndrho => dummy
     156        1224 :       e_rho_rho => dummy
     157        1224 :       e_ndrho_rho => dummy
     158        1224 :       e_ndrho_ndrho => dummy
     159        1224 :       e_rho_rho_rho => dummy
     160        1224 :       e_ndrho_rho_rho => dummy
     161        1224 :       e_ndrho_ndrho_rho => dummy
     162        1224 :       e_ndrho_ndrho_ndrho => dummy
     163             : 
     164        1224 :       IF (grad_deriv >= 0) THEN
     165             :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     166        1224 :                                          allocate_deriv=.TRUE.)
     167        1224 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     168             :       END IF
     169        1224 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     170             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     171        1216 :                                          allocate_deriv=.TRUE.)
     172        1216 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     173             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     174        1216 :                                          allocate_deriv=.TRUE.)
     175        1216 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     176             :       END IF
     177        1224 :       IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     178             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     179         378 :                                          allocate_deriv=.TRUE.)
     180         378 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     181             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
     182         378 :                                          allocate_deriv=.TRUE.)
     183         378 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
     184             :          deriv => xc_dset_get_derivative(deriv_set, &
     185         378 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     186         378 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     187             :       END IF
     188        1224 :       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        1224 :       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        1224 : !$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        1224 :       CALL timestop(handle)
     225        1224 :    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        1224 :    SUBROUTINE xb88_lr_lda_calc(rho, norm_drho, &
     255        1224 :                                e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
     256        1224 :                                e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
     257        1224 :                                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        1224 :       my_epsilon_rho = epsilon_rho
     291        1224 :       epsilon_rho43 = my_epsilon_rho**(4.0_dp/3.0_dp)
     292        1224 :       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    22177632 :          my_rho = rho(ii)*0.5_dp
     298    22177632 :          my_ndrho = norm_drho(ii)*0.5_dp
     299    22177632 :          IF (my_rho > my_epsilon_rho) THEN
     300    21158517 :             IF (grad_deriv >= 0) THEN
     301    21158517 :                t1 = my_rho**(0.1e1_dp/0.3e1_dp)
     302    21158517 :                t2 = t1*my_rho
     303    21158517 :                t3 = 0.1e1_dp/t2
     304    21158517 :                xx = my_ndrho*MAX(t3, epsilon_rho43)
     305    21158517 :                t5 = my_ndrho**2
     306    21158517 :                t6 = beta*t5
     307    21158517 :                t7 = my_rho**2
     308    21158517 :                t8 = t1**2
     309    21158517 :                t10 = 0.1e1_dp/t8/t7
     310    21158517 :                t11 = beta*my_ndrho
     311    21158517 :                t12 = LOG(xx + SQRT(xx**0.2e1_dp + 0.1e1_dp))
     312    21158517 :                t16 = 0.10e1_dp + 0.60e1_dp*t11*t3*t12
     313    21158517 :                t17 = 0.1e1_dp/t16
     314    21158517 :                t18 = t10*t17
     315    21158517 :                t21 = 0.20e1_dp*Cx + 0.20e1_dp*t6*t18
     316    21158517 :                t22 = SQRT(t21)
     317    21158517 :                t23 = t22*t21
     318    21158517 :                t24 = my_rho*t23
     319    21158517 :                t26 = rootpi
     320    21158517 :                t27 = 0.1e1_dp/t26
     321    21158517 :                t28 = 0.1e1_dp/omega
     322    21158517 :                t29 = 0.1e1_dp/t22
     323    21158517 :                t30 = t28*t29
     324    21158517 :                t31 = t26*t1
     325    21158517 :                t34 = erf(0.3000000000e1_dp*t30*t31)
     326    21158517 :                t36 = omega*t22
     327    21158517 :                t37 = 0.1e1_dp/t1
     328    21158517 :                t38 = t27*t37
     329    21158517 :                t39 = omega**2
     330    21158517 :                t40 = 0.1e1_dp/t39
     331    21158517 :                t41 = 0.1e1_dp/t21
     332    21158517 :                t42 = t40*t41
     333    21158517 :                t43 = pi*t8
     334    21158517 :                t44 = t42*t43
     335    21158517 :                t46 = EXP(-0.8999999998e1_dp*t44)
     336    21158517 :                t47 = t39*t21
     337    21158517 :                t48 = 0.1e1_dp/pi
     338    21158517 :                t49 = 0.1e1_dp/t8
     339    21158517 :                t51 = t46 - 0.10e1_dp
     340    21158517 :                t52 = t48*t49*t51
     341    21158517 :                t55 = t46 - 0.15e1_dp - 0.5555555558e-1_dp*t47*t52
     342    21158517 :                t59 = t26*t34 + 0.3333333334e0_dp*t36*t38*t55
     343    21158517 :                t60 = t27*t59
     344             : 
     345             :                !! Multiply with 2.0 because it code comes from LSD
     346    21158517 :                e_0(ii) = e_0(ii) - 0.2222222224e0_dp*t24*omega*t60*sx*2.0_dp
     347             : 
     348             :             END IF
     349    21158517 :             IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     350    20913217 :                t64 = t23*omega
     351    20913217 :                t68 = my_rho*t22*omega
     352    20913217 :                t69 = t7*my_rho
     353    20913217 :                t71 = 0.1e1_dp/t8/t69
     354    20913217 :                t72 = t71*t17
     355    20913217 :                t75 = t16**2
     356    20913217 :                t76 = 0.1e1_dp/t75
     357    20913217 :                t77 = t10*t76
     358    20913217 :                t79 = 0.1e1_dp/t1/t7
     359    20913217 :                t84 = 1 + t5*t10
     360    20913217 :                t85 = SQRT(t84)
     361    20913217 :                t86 = 0.1e1_dp/t85
     362    20913217 :                t87 = t71*t86
     363    20913217 :                t90 = -0.8000000000e1_dp*t11*t79*t12 - 0.8000000000e1_dp*t6*t87
     364    20913217 :                t91 = t77*t90
     365    20913217 :                t94 = -0.5333333333e1_dp*t6*t72 - 0.20e1_dp*t6*t91
     366    20913217 :                t95 = t60*t94
     367    20913217 :                t98 = omega*t27
     368    20913217 :                t99 = SQRT(0.3141592654e1_dp)
     369    20913217 :                t100 = 0.1e1_dp/t99
     370    20913217 :                t101 = t26*t100
     371    20913217 :                t103 = EXP(-0.9000000000e1_dp*t44)
     372    20913217 :                t104 = 0.1e1_dp/t23
     373    20913217 :                t105 = t28*t104
     374    20913217 :                t109 = t26*t49
     375             :                t112 = -0.1500000000e1_dp*t105*t31*t94 + 0.1000000000e1_dp*t30* &
     376    20913217 :                       t109
     377    20913217 :                t113 = t103*t112
     378    20913217 :                t116 = omega*t29
     379    20913217 :                t117 = t116*t27
     380    20913217 :                t118 = t37*t55
     381    20913217 :                t122 = t27*t3
     382    20913217 :                t126 = t21**2
     383    20913217 :                t127 = 0.1e1_dp/t126
     384    20913217 :                t128 = t40*t127
     385    20913217 :                t130 = t128*t43*t94
     386    20913217 :                t132 = pi*t37
     387    20913217 :                t133 = t42*t132
     388    20913217 :                t135 = 0.8999999998e1_dp*t130 - 0.5999999999e1_dp*t133
     389    20913217 :                t136 = t135*t46
     390    20913217 :                t137 = t39*t94
     391    20913217 :                t140 = t8*my_rho
     392    20913217 :                t141 = 0.1e1_dp/t140
     393    20913217 :                t143 = t48*t141*t51
     394    20913217 :                t146 = t47*t48
     395    20913217 :                t147 = t49*t135
     396    20913217 :                t148 = t147*t46
     397             :                t151 = t136 - 0.5555555558e-1_dp*t137*t52 + 0.3703703705e-1_dp*t47 &
     398    20913217 :                       *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    20913217 :                       t151
     402             : 
     403             :                e_rho(ii) = e_rho(ii) + (-0.2222222224e0_dp*t64*t60 - 0.3333333336e0_dp*t68*t95 - &
     404    20913217 :                                         0.2222222224e0_dp*t24*t98*t155)*sx
     405             : 
     406    20913217 :                t168 = 0.60e1_dp*beta*t3*t12 + 0.60e1_dp*t11*t10*t86
     407    20913217 :                t169 = t77*t168
     408    20913217 :                t172 = 0.40e1_dp*t11*t18 - 0.20e1_dp*t6*t169
     409    20913217 :                t173 = t60*t172
     410    20913217 :                t176 = pi*t100
     411    20913217 :                t177 = t176*t103
     412    20913217 :                t185 = t128*pi
     413    20913217 :                t186 = t8*t172
     414    20913217 :                t190 = t39*t172
     415             :                t196 = 0.8999999998e1_dp*t185*t186*t46 - 0.5555555558e-1_dp*t190 &
     416    20913217 :                       *t52 - 0.5000000001e0_dp*t41*t172*t46
     417             :                t200 = -0.3000000000e1_dp*t177*t105*t1*t172 + 0.1666666667e0_dp* &
     418    20913217 :                       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    20913217 :                                             *t200)*sx
     422             : 
     423             :             END IF
     424    21158517 :             IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     425     3723264 :                t207 = t27*t155
     426     3723264 :                t211 = my_rho*t29*omega
     427     3723264 :                t212 = t94**2
     428     3723264 :                t213 = t60*t212
     429     3723264 :                t216 = t207*t94
     430     3723264 :                t219 = t7**2
     431     3723264 :                t221 = 0.1e1_dp/t8/t219
     432     3723264 :                t222 = t221*t17
     433     3723264 :                t225 = t71*t76
     434     3723264 :                t226 = t225*t90
     435     3723264 :                t230 = 0.1e1_dp/t75/t16
     436     3723264 :                t231 = t10*t230
     437     3723264 :                t232 = t90**2
     438     3723264 :                t233 = t231*t232
     439     3723264 :                t237 = 0.1e1_dp/t1/t69
     440     3723264 :                t241 = t221*t86
     441     3723264 :                t244 = t5**2
     442     3723264 :                t245 = beta*t244
     443     3723264 :                t250 = 0.1e1_dp/t85/t84
     444     3723264 :                t251 = 0.1e1_dp/t1/t219/t69*t250
     445             :                t254 = 0.1866666667e2_dp*t11*t237*t12 + 0.4000000000e2_dp*t6*t241 &
     446     3723264 :                       - 0.1066666667e2_dp*t245*t251
     447     3723264 :                t255 = t77*t254
     448             :                t258 = 0.1955555555e2_dp*t6*t222 + 0.1066666667e2_dp*t6*t226 + 0.40e1_dp &
     449     3723264 :                       *t6*t233 - 0.20e1_dp*t6*t255
     450     3723264 :                t259 = t60*t258
     451     3723264 :                t264 = 0.9000000000e1_dp*t130 - 0.6000000000e1_dp*t133
     452     3723264 :                t265 = t264*t103
     453     3723264 :                t270 = 0.1e1_dp/t22/t126
     454     3723264 :                t271 = t28*t270
     455     3723264 :                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     3723264 :                       *t30*t281
     459     3723264 :                t285 = t103*t284
     460     3723264 :                t289 = omega*t104*t27
     461     3723264 :                t293 = t3*t55
     462     3723264 :                t297 = t37*t151
     463     3723264 :                t304 = t27*t79
     464     3723264 :                t311 = t126*t21
     465     3723264 :                t312 = 0.1e1_dp/t311
     466     3723264 :                t313 = t40*t312
     467     3723264 :                t316 = 0.1800000000e2_dp*t313*t43*t212
     468     3723264 :                t319 = 0.1200000000e2_dp*t128*t132*t94
     469     3723264 :                t321 = t128*t43*t258
     470     3723264 :                t323 = pi*t3
     471     3723264 :                t325 = 0.2000000000e1_dp*t42*t323
     472     3723264 :                t326 = -t316 + t319 + 0.8999999998e1_dp*t321 + t325
     473     3723264 :                t328 = t135**2
     474     3723264 :                t330 = t39*t258
     475     3723264 :                t335 = t137*t48
     476     3723264 :                t339 = t48*t10*t51
     477     3723264 :                t343 = t141*t135*t46
     478     3723264 :                t346 = t49*t326
     479     3723264 :                t347 = t346*t46
     480     3723264 :                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     3723264 :                       *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     3723264 :                       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     3723264 :                                                 *t68*t259 - 0.2222222224e0_dp*t24*t98*t358)*sx
     494             : 
     495     3723264 :                t365 = t27*t200
     496     3723264 :                t368 = t94*t172
     497     3723264 :                t372 = t365*t94
     498     3723264 :                t377 = t225*t168
     499     3723264 :                t382 = t6*t10
     500     3723264 :                t383 = t230*t90
     501     3723264 :                t384 = t383*t168
     502     3723264 :                t393 = beta*t5*my_ndrho
     503     3723264 :                t397 = 0.1e1_dp/t1/t219/t7*t250
     504             :                t400 = -0.8000000000e1_dp*beta*t79*t12 - 0.2400000000e2_dp*t11* &
     505     3723264 :                       t87 + 0.8000000000e1_dp*t393*t397
     506     3723264 :                t401 = t77*t400
     507             :                t404 = -0.1066666667e2_dp*t11*t72 + 0.5333333333e1_dp*t6*t377 - 0.40e1_dp &
     508     3723264 :                       *t11*t91 + 0.40e1_dp*t382*t384 - 0.20e1_dp*t6*t401
     509     3723264 :                t405 = t60*t404
     510     3723264 :                t408 = t207*t172
     511     3723264 :                t412 = t26*pi*t100
     512     3723264 :                t413 = t412*t128
     513     3723264 :                t417 = t271*t26
     514     3723264 :                t418 = t1*t94
     515             :                t428 = 0.2250000000e1_dp*t417*t418*t172 - 0.1500000000e1_dp*t105 &
     516     3723264 :                       *t31*t404 - 0.5000000000e0_dp*t105*t109*t172
     517     3723264 :                t429 = t103*t428
     518     3723264 :                t435 = t37*t196
     519     3723264 :                t451 = t313*pi
     520     3723264 :                t452 = t8*t94
     521     3723264 :                t455 = 0.1800000000e2_dp*t451*t452*t172
     522     3723264 :                t457 = t128*t43*t404
     523     3723264 :                t460 = t128*t132*t172
     524     3723264 :                t462 = -t455 + 0.8999999998e1_dp*t457 + 0.5999999999e1_dp*t460
     525     3723264 :                t463 = t462*t46
     526     3723264 :                t464 = t135*t40
     527     3723264 :                t465 = t464*t127
     528     3723264 :                t466 = t172*t46
     529     3723264 :                t467 = t43*t466
     530     3723264 :                t470 = t39*t404
     531     3723264 :                t473 = t94*t127
     532     3723264 :                t478 = 0.1e1_dp/my_rho
     533     3723264 :                t479 = t41*t478
     534     3723264 :                t482 = t190*t48
     535     3723264 :                t486 = t49*t462*t46
     536     3723264 :                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     3723264 :                       - 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     3723264 :                       *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     3723264 :                                                     *t24*t98*t496)*sx
     551             : 
     552     3723264 :                t501 = t172**2
     553     3723264 :                t502 = t60*t501
     554     3723264 :                t505 = t365*t172
     555     3723264 :                t508 = beta*t10
     556     3723264 :                t513 = t168**2
     557     3723264 :                t514 = t231*t513
     558     3723264 :                t519 = t219*my_rho
     559     3723264 :                t521 = 0.1e1_dp/t1/t519
     560     3723264 :                t525 = 0.120e2_dp*t508*t86 - 0.60e1_dp*t6*t521*t250
     561     3723264 :                t526 = t77*t525
     562             :                t529 = 0.40e1_dp*t508*t17 - 0.80e1_dp*t11*t169 + 0.40e1_dp*t6*t514 &
     563     3723264 :                       - 0.20e1_dp*t6*t526
     564     3723264 :                t530 = t60*t529
     565     3723264 :                t533 = pi**2
     566     3723264 :                t534 = t533*t100
     567     3723264 :                t536 = 0.1e1_dp/t39/omega
     568     3723264 :                t537 = t534*t536
     569     3723264 :                t539 = 0.1e1_dp/t22/t311
     570     3723264 :                t562 = t8*t501
     571     3723264 :                t566 = t8*t529
     572     3723264 :                t570 = t39**2
     573     3723264 :                t571 = 0.1e1_dp/t570
     574     3723264 :                t572 = t126**2
     575     3723264 :                t573 = 0.1e1_dp/t572
     576     3723264 :                t574 = t571*t573
     577     3723264 :                t575 = t574*t533
     578     3723264 :                t576 = t2*t501
     579     3723264 :                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     3723264 :                       *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     3723264 :                       *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     3723264 :                                    *sx
     592             : 
     593             :             END IF
     594    21158517 :             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        1224 :    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 1.15