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

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief calculates the pbe correlation functional
      10              : !> \note
      11              : !>      This was generated with the help of a maple worksheet that you can
      12              : !>      find under doc/pbe.mw .
      13              : !>      I did not add 3. derivatives in the polazied (lsd) case because it
      14              : !>      would have added 2500 lines of code. If you need them the worksheet
      15              : !>      is already prepared for them, and by uncommenting a couple of lines
      16              : !>      you should be able to generate them.
      17              : !> \par History
      18              : !>      09.2004 created [fawzi]
      19              : !> \author fawzi
      20              : ! **************************************************************************************************
      21              : MODULE xc_pbe
      22              :    USE bibliography,                    ONLY: Perdew1996,&
      23              :                                               Perdew2008,&
      24              :                                               Zhang1998,&
      25              :                                               cite_reference
      26              :    USE input_section_types,             ONLY: section_vals_type,&
      27              :                                               section_vals_val_get
      28              :    USE kinds,                           ONLY: dp
      29              :    USE mathconstants,                   ONLY: pi
      30              :    USE xc_derivative_desc,              ONLY: deriv_norm_drho,&
      31              :                                               deriv_norm_drhoa,&
      32              :                                               deriv_norm_drhob,&
      33              :                                               deriv_rho,&
      34              :                                               deriv_rhoa,&
      35              :                                               deriv_rhob
      36              :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      37              :                                               xc_dset_get_derivative
      38              :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      39              :                                               xc_derivative_type
      40              :    USE xc_input_constants,              ONLY: xc_pbe_orig,&
      41              :                                               xc_pbe_rev,&
      42              :                                               xc_pbe_sol
      43              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      44              :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      45              :                                               xc_rho_set_type
      46              : #include "../base/base_uses.f90"
      47              : 
      48              :    IMPLICIT NONE
      49              :    PRIVATE
      50              : 
      51              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      52              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pbe'
      53              :    REAL(kind=dp), PARAMETER, PRIVATE :: a = 0.04918_dp, b = 0.132_dp, &
      54              :                                         c = 0.2533_dp, d = 0.349_dp
      55              : 
      56              :    PUBLIC :: pbe_lda_info, pbe_lsd_info, pbe_lda_eval, pbe_lsd_eval
      57              : CONTAINS
      58              : 
      59              : ! **************************************************************************************************
      60              : !> \brief return various information on the functional
      61              : !> \param pbe_params section selecting the various parameters for the functional
      62              : !> \param reference string with the reference of the actual functional
      63              : !> \param shortform string with the shortform of the functional name
      64              : !> \param needs the components needed by this functional are set to
      65              : !>        true (does not set the unneeded components to false)
      66              : !> \param max_deriv ...
      67              : !> \author fawzi
      68              : ! **************************************************************************************************
      69       255228 :    SUBROUTINE pbe_lda_info(pbe_params, reference, shortform, needs, max_deriv)
      70              :       TYPE(section_vals_type), POINTER                   :: pbe_params
      71              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      72              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      73              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      74              : 
      75              :       INTEGER                                            :: param
      76              :       REAL(kind=dp)                                      :: sc, sx
      77              : 
      78        85076 :       CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
      79        85076 :       CALL section_vals_val_get(pbe_params, "scale_x", r_val=sx)
      80        85076 :       CALL section_vals_val_get(pbe_params, "scale_c", r_val=sc)
      81              : 
      82        84723 :       SELECT CASE (param)
      83              :       CASE (xc_pbe_orig)
      84        84723 :          CALL cite_reference(Perdew1996)
      85        84723 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
      86        55185 :             IF (PRESENT(reference)) THEN
      87              :                reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "// &
      88          409 :                            "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996) {spin unpolarized}"
      89              :             END IF
      90        55185 :             IF (PRESENT(shortform)) THEN
      91          409 :                shortform = "PBE Perdew-Burke-Ernzerhof xc functional (unpolarized)"
      92              :             END IF
      93              :          ELSE
      94        29538 :             IF (PRESENT(reference)) THEN
      95              :                WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
      96          169 :                   "J.P.Perdew, K.Burke, M.Ernzerhof, ", &
      97          169 :                   "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)", &
      98          338 :                   sx, sc, "{spin unpolarized}"
      99              :             END IF
     100        29538 :             IF (PRESENT(shortform)) THEN
     101              :                WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") &
     102          169 :                   "PBE, Perdew-Burke-Ernzerhof xc functional (unpolarized)", sx, sc
     103              :             END IF
     104              :          END IF
     105              :       CASE (xc_pbe_rev)
     106          157 :          CALL cite_reference(Perdew1996)
     107          157 :          CALL cite_reference(Zhang1998)
     108          157 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
     109            0 :             IF (PRESENT(reference)) THEN
     110              :                reference = "revPBE, Yingkay Zhang and Weitao Yang,"// &
     111            0 :                            " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998) {spin unpolarized}"
     112              :             END IF
     113            0 :             IF (PRESENT(shortform)) THEN
     114            0 :                shortform = "revPBE, revised PBE xc functional (unpolarized)"
     115              :             END IF
     116              :          ELSE
     117          157 :             IF (PRESENT(reference)) THEN
     118              :                WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
     119            3 :                   "revPBE, Yingkay Zhang and Weitao Yang,", &
     120            3 :                   " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)", &
     121            6 :                   sx, sc, "{spin unpolarized}"
     122              :             END IF
     123          157 :             IF (PRESENT(shortform)) THEN
     124              :                WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") &
     125            3 :                   "revPBE, revised PBE xc functional (unpolarized)", sx, sc
     126              :             END IF
     127              :          END IF
     128              :       CASE (xc_pbe_sol)
     129          196 :          CALL cite_reference(Perdew1996)
     130          196 :          CALL cite_reference(Perdew2008)
     131          196 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
     132          166 :             IF (PRESENT(reference)) THEN
     133              :                reference = "PBEsol, J.P. Perdew et al., "// &
     134              :                            "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) "// &
     135            2 :                            "{spin unpolarized}"
     136              :             END IF
     137          166 :             IF (PRESENT(shortform)) THEN
     138            2 :                shortform = "PBEsol, PBE for solids and surfaces xc functional (unpolarized)"
     139              :             END IF
     140              :          ELSE
     141           30 :             IF (PRESENT(reference)) THEN
     142              :                WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
     143            2 :                   "PBEsol, J.P. Perdew et al., ", &
     144            2 :                   "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) ", &
     145            4 :                   sx, sc, "{spin unpolarized}"
     146              :             END IF
     147           30 :             IF (PRESENT(shortform)) THEN
     148              :                WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") &
     149            2 :                   "PBEsol, PBE for solids and surfaces xc functional (unpolarized)", sx, sc
     150              :             END IF
     151              :          END IF
     152              :       CASE default
     153        85076 :          CPABORT("Unsupported parametrization")
     154              :       END SELECT
     155        85076 :       IF (PRESENT(needs)) THEN
     156        84491 :          needs%rho = .TRUE.
     157        84491 :          needs%norm_drho = .TRUE.
     158              :       END IF
     159        85076 :       IF (PRESENT(max_deriv)) max_deriv = 3
     160              : 
     161        85076 :    END SUBROUTINE pbe_lda_info
     162              : 
     163              : ! **************************************************************************************************
     164              : !> \brief return various information on the functional
     165              : !> \param pbe_params section selecting the various parameters for the functional
     166              : !> \param reference string with the reference of the actual functional
     167              : !> \param shortform string with the shortform of the functional name
     168              : !> \param needs the components needed by this functional are set to
     169              : !>        true (does not set the unneeded components to false)
     170              : !> \param max_deriv ...
     171              : !> \author fawzi
     172              : ! **************************************************************************************************
     173        57669 :    SUBROUTINE pbe_lsd_info(pbe_params, reference, shortform, needs, max_deriv)
     174              :       TYPE(section_vals_type), POINTER                   :: pbe_params
     175              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
     176              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
     177              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
     178              : 
     179              :       INTEGER                                            :: param
     180              :       REAL(kind=dp)                                      :: sc, sx
     181              : 
     182        19223 :       CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
     183        19223 :       CALL section_vals_val_get(pbe_params, "scale_x", r_val=sx)
     184        19223 :       CALL section_vals_val_get(pbe_params, "scale_c", r_val=sc)
     185              : 
     186        19051 :       SELECT CASE (param)
     187              :       CASE (xc_pbe_orig)
     188        19051 :          CALL cite_reference(Perdew1996)
     189        19051 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
     190        13459 :             IF (PRESENT(reference)) THEN
     191              :                reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "// &
     192          248 :                            "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996) {spin polarized}"
     193              :             END IF
     194        13459 :             IF (PRESENT(shortform)) THEN
     195          248 :                shortform = "PBE xc functional (spin polarized)"
     196              :             END IF
     197              :          ELSE
     198         5592 :             IF (PRESENT(reference)) THEN
     199              :                WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
     200           50 :                   "J.P.Perdew, K.Burke, M.Ernzerhof, ", &
     201           50 :                   "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)", &
     202          100 :                   sx, sc, "{spin polarized}"
     203              :             END IF
     204         5592 :             IF (PRESENT(shortform)) THEN
     205              :                WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") &
     206           50 :                   "PBE xc functional (spin polarized)", sx, sc
     207              :             END IF
     208              :          END IF
     209              :       CASE (xc_pbe_rev)
     210           24 :          CALL cite_reference(Perdew1996)
     211           24 :          CALL cite_reference(Zhang1998)
     212           24 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
     213            0 :             IF (PRESENT(reference)) THEN
     214              :                reference = "revPBE, Yingkay Zhang and Weitao Yang,"// &
     215            0 :                            " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998) {spin polarized}"
     216              :             END IF
     217            0 :             IF (PRESENT(shortform)) THEN
     218            0 :                shortform = "revPBE, revised PBE xc functional (spin polarized)"
     219              :             END IF
     220              :          ELSE
     221           24 :             IF (PRESENT(reference)) THEN
     222              :                WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
     223            0 :                   "revPBE, Yingkay Zhang and Weitao Yang,", &
     224            0 :                   " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)", &
     225            0 :                   sx, sc, "{spin polarized}"
     226              :             END IF
     227           24 :             IF (PRESENT(shortform)) THEN
     228              :                WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") &
     229            0 :                   "revPBE, revised PBE xc functional (spin polarized)", sx, sc
     230              :             END IF
     231              :          END IF
     232              :       CASE (xc_pbe_sol)
     233          148 :          CALL cite_reference(Perdew1996)
     234          148 :          CALL cite_reference(Perdew2008)
     235          148 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
     236          145 :             IF (PRESENT(reference)) THEN
     237              :                reference = "PBEsol, J.P. Perdew et al., "// &
     238              :                            "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) "// &
     239            1 :                            "{spin polarized}"
     240              :             END IF
     241          145 :             IF (PRESENT(shortform)) THEN
     242            1 :                shortform = "PBEsol, PBE for solids and surfaces xc functional (spin polarized)"
     243              :             END IF
     244              :          ELSE
     245            3 :             IF (PRESENT(reference)) THEN
     246              :                WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
     247            1 :                   "PBEsol, J.P. Perdew et al., ", &
     248            1 :                   "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) ", &
     249            2 :                   sx, sc, "{spin unpolarized}"
     250              :             END IF
     251            3 :             IF (PRESENT(shortform)) THEN
     252              :                WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") &
     253            1 :                   "PBEsol, PBE for solids and surfaces xc functional (spin polarized)", sx, sc
     254              :             END IF
     255              :          END IF
     256              :       CASE default
     257        19223 :          CPABORT("Unsupported parametrization")
     258              :       END SELECT
     259        19223 :       IF (PRESENT(needs)) THEN
     260        18923 :          needs%rho_spin = .TRUE.
     261        18923 :          needs%norm_drho_spin = .TRUE.
     262        18923 :          needs%norm_drho = .TRUE.
     263              :       END IF
     264        19223 :       IF (PRESENT(max_deriv)) max_deriv = 2
     265              : 
     266        19223 :    END SUBROUTINE pbe_lsd_info
     267              : 
     268              : ! **************************************************************************************************
     269              : !> \brief evaluates the pbe correlation functional for lda
     270              : !> \param rho_set the density where you want to evaluate the functional
     271              : !> \param deriv_set place where to store the functional derivatives (they are
     272              : !>        added to the derivatives)
     273              : !> \param grad_deriv degree of the derivative that should be evaluated,
     274              : !>        if positive all the derivatives up to the given degree are evaluated,
     275              : !>        if negative only the given degree is calculated
     276              : !> \param pbe_params ...
     277              : !> \author fawzi
     278              : ! **************************************************************************************************
     279       448495 :    SUBROUTINE pbe_lda_eval(rho_set, deriv_set, grad_deriv, pbe_params)
     280              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     281              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     282              :       INTEGER, INTENT(in)                                :: grad_deriv
     283              :       TYPE(section_vals_type), POINTER                   :: pbe_params
     284              : 
     285              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pbe_lda_eval'
     286              : 
     287              :       INTEGER                                            :: handle, npoints, param
     288              :       INTEGER, DIMENSION(2, 3)                           :: bo
     289              :       REAL(kind=dp)                                      :: epsilon_rho, scale_ec, scale_ex
     290        89699 :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
     291        89699 :          e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
     292        89699 :          e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
     293              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     294              : 
     295        89699 :       CALL timeset(routineN, handle)
     296              : 
     297              :       CALL xc_rho_set_get(rho_set, rho=rho, &
     298        89699 :                           norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
     299        89699 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     300              : 
     301        89699 :       dummy => rho
     302              : 
     303        89699 :       e_0 => dummy
     304        89699 :       e_rho => dummy
     305        89699 :       e_ndrho => dummy
     306        89699 :       e_rho_rho => dummy
     307        89699 :       e_ndrho_rho => dummy
     308        89699 :       e_ndrho_ndrho => dummy
     309        89699 :       e_rho_rho_rho => dummy
     310        89699 :       e_ndrho_rho_rho => dummy
     311        89699 :       e_ndrho_ndrho_rho => dummy
     312        89699 :       e_ndrho_ndrho_ndrho => dummy
     313              : 
     314        89699 :       IF (grad_deriv >= 0) THEN
     315              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     316        89699 :                                          allocate_deriv=.TRUE.)
     317        89699 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     318              :       END IF
     319        89699 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     320              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     321        87167 :                                          allocate_deriv=.TRUE.)
     322        87167 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     323              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     324        87167 :                                          allocate_deriv=.TRUE.)
     325        87167 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     326              :       END IF
     327        89699 :       IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     328              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     329        11724 :                                          allocate_deriv=.TRUE.)
     330        11724 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     331              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
     332        11724 :                                          allocate_deriv=.TRUE.)
     333        11724 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
     334              :          deriv => xc_dset_get_derivative(deriv_set, &
     335        11724 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     336        11724 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     337              :       END IF
     338        89699 :       IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
     339              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
     340            0 :                                          allocate_deriv=.TRUE.)
     341            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     342              :          deriv => xc_dset_get_derivative(deriv_set, &
     343            0 :                                          [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.TRUE.)
     344            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
     345              :          deriv => xc_dset_get_derivative(deriv_set, &
     346            0 :                                          [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.TRUE.)
     347            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
     348              :          deriv => xc_dset_get_derivative(deriv_set, &
     349            0 :                                          [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     350            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
     351              :       END IF
     352        89699 :       IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
     353            0 :          CPABORT("derivatives bigger than 3 not implemented")
     354              :       END IF
     355              : 
     356        89699 :       CALL section_vals_val_get(pbe_params, "scale_c", r_val=scale_ec)
     357        89699 :       CALL section_vals_val_get(pbe_params, "scale_x", r_val=scale_ex)
     358        89699 :       CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
     359              : 
     360              : !$OMP PARALLEL DEFAULT(NONE), &
     361              : !$OMP SHARED(rho,norm_drho,e_0,e_rho,e_ndrho,e_rho_rho,e_ndrho_rho),&
     362              : !$OMP SHARED(e_ndrho_ndrho,e_rho_rho_rho,e_ndrho_rho_rho,e_ndrho_ndrho_rho),&
     363              : !$OMP SHARED(e_ndrho_ndrho_ndrho,grad_deriv,npoints,epsilon_rho),&
     364              : !$OMP SHARED(pbe_params),&
     365        89699 : !$OMP SHARED(param,scale_ec,scale_ex)
     366              :       CALL pbe_lda_calc(rho=rho, norm_drho=norm_drho, &
     367              :                         e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
     368              :                         e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
     369              :                         e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
     370              :                         e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, &
     371              :                         grad_deriv=grad_deriv, &
     372              :                         npoints=npoints, epsilon_rho=epsilon_rho, &
     373              :                         param=param, scale_ec=scale_ec, scale_ex=scale_ex)
     374              : !$OMP END PARALLEL
     375              : 
     376        89699 :       CALL timestop(handle)
     377        89699 :    END SUBROUTINE pbe_lda_eval
     378              : 
     379              : ! **************************************************************************************************
     380              : !> \brief evaluates the pbe correlation functional for lda
     381              : !> \param rho the density where you want to evaluate the functional
     382              : !> \param norm_drho ...
     383              : !> \param e_0 ...
     384              : !> \param e_rho ...
     385              : !> \param e_ndrho ...
     386              : !> \param e_rho_rho ...
     387              : !> \param e_ndrho_rho ...
     388              : !> \param e_ndrho_ndrho ...
     389              : !> \param e_rho_rho_rho ...
     390              : !> \param e_ndrho_rho_rho ...
     391              : !> \param e_ndrho_ndrho_rho ...
     392              : !> \param e_ndrho_ndrho_ndrho ...
     393              : !> \param grad_deriv degree of the derivative that should be evaluated,
     394              : !>        if positive all the derivatives up to the given degree are evaluated,
     395              : !>        if negative only the given degree is calculated
     396              : !> \param npoints ...
     397              : !> \param epsilon_rho ...
     398              : !> \param ! ...
     399              : !> \param pbe_params parameters for the pbe functional
     400              : !> \author fawzi
     401              : ! **************************************************************************************************
     402        89699 :    SUBROUTINE pbe_lda_calc(rho, norm_drho, &
     403        89699 :                            e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
     404        89699 :                            e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
     405        89699 :                            e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, &
     406              :                            !     pbe_params)
     407              :                            param, scale_ec, scale_ex)
     408              :       INTEGER, INTENT(in)                                :: npoints, grad_deriv
     409              :       REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_ndrho, &
     410              :          e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, &
     411              :          e_ndrho, e_rho, e_0
     412              :       REAL(kind=dp), DIMENSION(1:npoints), INTENT(in)    :: norm_drho, rho
     413              :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho
     414              :       INTEGER, INTENT(in)                                :: param
     415              :       REAL(kind=dp), INTENT(in)                          :: scale_ec, scale_ex
     416              : 
     417              :       INTEGER                                            :: ii
     418              :       REAL(kind=dp) :: A, A1rho, A1rhorho, A2rho, A_1, alpha_1_1, Arho, Arho1rho, Arhorho, &
     419              :          Arhorhorho, beta, beta_1_1, beta_2_1, beta_3_1, beta_4_1, e_c_u_0, e_c_u_01rho, &
     420              :          e_c_u_01rhorho, e_c_u_02rho, e_c_u_0rho, e_c_u_0rho1rho, e_c_u_0rhorho, e_c_u_0rhorhorho, &
     421              :          epsilon_cGGA, epsilon_cGGArho, epsilon_cGGArhorho, ex_ldarhorhorho, ex_unif, ex_unif1rho, &
     422              :          ex_unif1rhorho, ex_unif2rho, ex_unifrho, ex_unifrho1rho, ex_unifrhorho, Fx, Fx1rho, &
     423              :          Fx1rhorho, Fx2rho, Fxnorm_drho, Fxnorm_drho1rho, Fxnorm_drhonorm_drho, Fxnorm_drhorho, &
     424              :          Fxrho, Fxrho1rho, Fxrhorho, gamma_var, Hnorm_drho, Hnorm_drhonorm_drho
     425              :       REAL(kind=dp) :: Hnorm_drhorho, k_f, k_f2rho, k_frho, k_frhorho, k_frhorhorho, k_s, k_s1rho, &
     426              :          k_s1rhorho, k_s2rho, k_srho, k_srho1rho, k_srhorho, kappa, kf, kf2rho, kfrho, kfrhorho, &
     427              :          kfrhorhorho, mu, my_norm_drho, my_rho, p_1, p_2, rs, rs2rho, rsrho, rsrhorho, &
     428              :          rsrhorhorho, s, s1rho, s1rhorho, s2norm_drho, s2rho, snorm_drho, snorm_drho1rho, &
     429              :          snorm_drhorho, srho, srho1rho, srhorho, t, t1, t1001, t1004, t1005, t1006, t1008, t101, &
     430              :          t1011, t1012, t1013, t1014, t1016, t1017, t1018, t1019, t1022, t1024, t1026, t1028, t103, &
     431              :          t1030, t1031, t1035, t1042, t1046, t1048, t105, t1050, t1052, t1054, t1055
     432              :       REAL(kind=dp) :: t1060, t1061, t1062, t1067, t108, t1093, t1097, t11, t1103, t1104, t1106, &
     433              :          t1109, t111, t1115, t1118, t1121, t1122, t1126, t1129, t113, t114, t1148, t115, t1152, &
     434              :          t1157, t1167, t1187, t119, t1196, t1203, t1218, t1226, t123, t124, t1249, t125, t126, &
     435              :          t1262, t1263, t127, t1291, t1292, t1295, t13, t131, t1329, t1342, t135, t138, t1380, &
     436              :          t1385, t1386, t1387, t1389, t139, t1393, t14, t140, t142, t146, t1468, t148, t1482, t149, &
     437              :          t1492, t1493, t150, t1504, t1505, t151, t1511, t1515, t1521, t1525, t1528, t153, t1532, &
     438              :          t1545, t155, t1565, t157, t1573, t158, t1584, t159, t1608, t1612
     439              :       REAL(kind=dp) :: t162, t1628, t1629, t163, t1633, t164, t1646, t1652, t1660, t167, t1672, &
     440              :          t1676, t168, t17, t170, t171, t1715, t1722, t175, t1758, t1759, t176, t1761, t177, t178, &
     441              :          t179, t1797, t182, t183, t1838, t1851, t186, t187, t1878, t1889, t189, t19, t190, t1907, &
     442              :          t191, t1922, t1927, t1933, t1937, t195, t1952, t196, t1964, t1965, t1968, t1969, t197, &
     443              :          t1972, t198, t1990, t1rho, t1rhorho, t2, t20, t200, t202, t2020, t2024, t2028, t2031, &
     444              :          t204, t2041, t208, t21, t214, t217, t218, t22, t226, t229, t230, t238, t239, t240, t241, &
     445              :          t246, t252, t253, t255, t259, t26, t260, t261, t262, t265, t266
     446              :       REAL(kind=dp) :: t267, t27, t271, t272, t273, t277, t278, t280, t281, t286, t290, t291, &
     447              :          t293, t294, t295, t296, t297, t299, t2norm_drho, t2rho, t3, t305, t309, t310, t315, t317, &
     448              :          t318, t321, t323, t324, t327, t329, t330, t331, t336, t338, t339, t340, t341, t348, t349, &
     449              :          t354, t357, t359, t360, t361, t362, t369, t370, t371, t374, t377, t378, t381, t382, t384, &
     450              :          t385, t387, t388, t390, t392, t393, t397, t4, t400, t401, t404, t408, t409, t410, t414, &
     451              :          t417, t418, t423, t426, t427, t432, t435, t436, t438, t439, t440, t448, t449, t451, t456, &
     452              :          t457, t458, t459, t461, t462, t463, t465, t466, t469, t470
     453              :       REAL(kind=dp) :: t471, t472, t476, t487, t491, t496, t498, t5, t500, t503, t506, t507, t510, &
     454              :          t513, t517, t521, t525, t529, t530, t535, t541, t548, t553, t556, t557, t559, t562, t566, &
     455              :          t581, t586, t590, t591, t594, t598, t6, t605, t609, t611, t612, t614, t627, t645, t65, &
     456              :          t654, t656, t661, t664, t669, t670, t671, t673, t675, t685, t69, t693, t7, t70, t701, &
     457              :          t702, t71, t714, t717, t72, t720, t73, t74, t740, t743, t748, t75, t758, t77, t776, t78, &
     458              :          t8, t80, t801, t809, t81, t812, t821, t825, t83, t831, t84, t85, t86, t868, t87, t877, &
     459              :          t879, t88, t880, t885, t90, t91, t94, t940, t942, t943, t945
     460              :       REAL(kind=dp) :: t946, t948, t95, t950, t951, t954, t967, t976, t98, t980, t982, t984, t989, &
     461              :          t99, t990, t994, t995, t998, t999, tnorm_drho, tnorm_drho1rho, tnorm_drhorho, &
     462              :          tnorm_drhorhorho, trho, trho1rho, trhorho, trhorhorho
     463              : 
     464              : !TYPE(section_vals_type), POINTER         :: pbe_params
     465              : !, param
     466              : ! scale_ec, scale_ex, snorm_drho, snorm_drho1rho, snorm_drhorho, srho, &
     467              : ! Parametrization of PBE
     468              : 
     469        89699 :       SELECT CASE (param)
     470              :          ! Original PBE
     471              :       CASE (xc_pbe_orig)
     472              :          kappa = 0.804e0_dp
     473              :          beta = 0.66725e-1_dp
     474          130 :          mu = beta*pi**2/0.3e1_dp
     475              :          ! Revised PBE (revPBE)
     476              :       CASE (xc_pbe_rev)
     477          130 :          kappa = 0.1245e1_dp
     478          130 :          beta = 0.66725e-1_dp
     479          130 :          mu = beta*pi**2/0.3e1_dp
     480              :          ! PBE for solids (PBEsol)
     481              :       CASE (xc_pbe_sol)
     482          188 :          kappa = 0.804e0_dp
     483          188 :          beta = 0.46e-1_dp
     484          188 :          mu = 0.1e2_dp/0.81e2_dp
     485              :       CASE default
     486        89699 :          CPABORT("Unsupported parametrization")
     487              :       END SELECT
     488              : 
     489        89699 :       gamma_var = (0.1e1_dp - LOG(0.2e1_dp))/pi**2
     490        89699 :       p_1 = 0.10e1_dp
     491        89699 :       A_1 = 0.31091e-1_dp
     492        89699 :       alpha_1_1 = 0.21370e0_dp
     493        89699 :       beta_1_1 = 0.75957e1_dp
     494        89699 :       beta_2_1 = 0.35876e1_dp
     495        89699 :       beta_3_1 = 0.16382e1_dp
     496        89699 :       beta_4_1 = 0.49294e0_dp
     497        89699 :       p_2 = 0.10e1_dp
     498        89699 :       t1 = 3**(0.1e1_dp/0.3e1_dp)
     499        89699 :       t2 = 4**(0.1e1_dp/0.3e1_dp)
     500        89699 :       t3 = t2**2
     501        89699 :       t4 = t1*t3
     502        89699 :       t5 = 0.1e1_dp/pi
     503        89699 :       t69 = pi**2
     504        89699 :       t627 = 0.1e1_dp/t69/pi
     505              : 
     506              :       SELECT CASE (grad_deriv)
     507              :       CASE default
     508              : !$OMP DO
     509              :          DO ii = 1, npoints
     510   1981970634 :             my_rho = rho(ii)
     511   1981970634 :             IF (my_rho > epsilon_rho) THEN
     512   1831948030 :                my_norm_drho = norm_drho(ii)
     513              : 
     514   1831948030 :                t6 = 0.1e1_dp/my_rho
     515   1831948030 :                t7 = t5*t6
     516   1831948030 :                t8 = t7**(0.1e1_dp/0.3e1_dp)
     517   1831948030 :                rs = t4*t8/0.4e1_dp
     518   1831948030 :                t11 = 0.1e1_dp + alpha_1_1*rs
     519   1831948030 :                t13 = 0.1e1_dp/A_1
     520   1831948030 :                t14 = SQRT(rs)
     521   1831948030 :                t17 = t14*rs
     522   1831948030 :                t19 = p_1 + 0.1e1_dp
     523   1831948030 :                t20 = rs**t19
     524   1831948030 :                t21 = beta_4_1*t20
     525   1831948030 :                t22 = beta_1_1*t14 + beta_2_1*rs + beta_3_1*t17 + t21
     526   1831948030 :                t26 = 0.1e1_dp + t13/t22/0.2e1_dp
     527   1831948030 :                t27 = LOG(t26)
     528   1831948030 :                e_c_u_0 = -0.2e1_dp*A_1*t11*t27
     529   1831948030 :                t65 = 2**(0.1e1_dp/0.3e1_dp)
     530   1831948030 :                t70 = t69*my_rho
     531   1831948030 :                t71 = t70**(0.1e1_dp/0.3e1_dp)
     532   1831948030 :                k_f = t1*t71
     533   1831948030 :                t72 = k_f*t5
     534   1831948030 :                t73 = SQRT(t72)
     535   1831948030 :                k_s = 0.2e1_dp*t73
     536   1831948030 :                t74 = 0.1e1_dp/k_s
     537   1831948030 :                t75 = my_norm_drho*t74
     538   1831948030 :                t = t75*t6/0.2e1_dp
     539   1831948030 :                t77 = 0.1e1_dp/gamma_var
     540   1831948030 :                t78 = beta*t77
     541   1831948030 :                t80 = EXP(-e_c_u_0*t77)
     542   1831948030 :                t81 = -0.1e1_dp + t80
     543   1831948030 :                A = t78/t81
     544   1831948030 :                t83 = t**2
     545   1831948030 :                t84 = A*t83
     546   1831948030 :                t85 = 0.1e1_dp + t84
     547   1831948030 :                t86 = t83*t85
     548   1831948030 :                t87 = A**2
     549   1831948030 :                t88 = t83**2
     550   1831948030 :                t90 = 0.1e1_dp + t84 + t87*t88
     551   1831948030 :                t91 = 0.1e1_dp/t90
     552   1831948030 :                t94 = 0.1e1_dp + t78*t86*t91
     553   1831948030 :                t95 = LOG(t94)
     554   1831948030 :                epsilon_cGGA = e_c_u_0 + gamma_var*t95
     555   1831948030 :                kf = k_f
     556   1831948030 :                ex_unif = -0.3e1_dp/0.4e1_dp*t5*kf
     557   1831948030 :                t98 = 0.1e1_dp/kf
     558   1831948030 :                t99 = my_norm_drho*t98
     559   1831948030 :                s = t99*t6/0.2e1_dp
     560   1831948030 :                t101 = s**2
     561   1831948030 :                t103 = 0.1e1_dp/kappa
     562   1831948030 :                t105 = 0.1e1_dp + mu*t101*t103
     563   1831948030 :                Fx = 0.1e1_dp + kappa - kappa/t105
     564   1831948030 :                t108 = my_rho*ex_unif
     565              : 
     566   1831948030 :                IF (grad_deriv >= 0) THEN
     567              :                   e_0(ii) = e_0(ii) + &
     568   1831948030 :                             scale_ex*t108*Fx + scale_ec*my_rho*epsilon_cGGA
     569              :                END IF
     570              : 
     571   1831948030 :                t111 = t8**2
     572   1831948030 :                t113 = 0.1e1_dp/t111*t5
     573   1831948030 :                t114 = my_rho**2
     574   1831948030 :                t115 = 0.1e1_dp/t114
     575   1831948030 :                rsrho = -t4*t113*t115/0.12e2_dp
     576   1831948030 :                t119 = A_1*alpha_1_1
     577   1831948030 :                t123 = t22**2
     578   1831948030 :                t124 = 0.1e1_dp/t123
     579   1831948030 :                t125 = t11*t124
     580   1831948030 :                t126 = 0.1e1_dp/t14
     581   1831948030 :                t127 = beta_1_1*t126
     582   1831948030 :                t131 = beta_3_1*t14
     583   1831948030 :                t135 = 0.1e1_dp/rs
     584              :                t138 = t127*rsrho/0.2e1_dp + beta_2_1*rsrho + 0.3e1_dp/ &
     585   1831948030 :                       0.2e1_dp*t131*rsrho + t21*t19*rsrho*t135
     586   1831948030 :                t139 = 0.1e1_dp/t26
     587   1831948030 :                t140 = t138*t139
     588   1831948030 :                e_c_u_0rho = -0.2e1_dp*t119*rsrho*t27 + t125*t140
     589   1831948030 :                t142 = t71**2
     590   1831948030 :                k_frho = t1/t142*t69/0.3e1_dp
     591   1831948030 :                t146 = 0.1e1_dp/t73
     592   1831948030 :                k_srho = t146*k_frho*t5
     593   1831948030 :                t148 = k_s**2
     594   1831948030 :                t149 = 0.1e1_dp/t148
     595   1831948030 :                t150 = my_norm_drho*t149
     596   1831948030 :                t151 = t6*k_srho
     597   1831948030 :                t153 = t75*t115
     598   1831948030 :                trho = -t150*t151/0.2e1_dp - t153/0.2e1_dp
     599   1831948030 :                t155 = gamma_var**2
     600   1831948030 :                t157 = beta/t155
     601   1831948030 :                t158 = t81**2
     602   1831948030 :                t159 = 0.1e1_dp/t158
     603   1831948030 :                Arho = t157*t159*e_c_u_0rho*t80
     604   1831948030 :                t162 = t78*t
     605   1831948030 :                t163 = t85*t91
     606   1831948030 :                t164 = t163*trho
     607   1831948030 :                t167 = Arho*t83
     608   1831948030 :                t168 = A*t
     609   1831948030 :                t170 = 0.2e1_dp*t168*trho
     610   1831948030 :                t171 = t167 + t170
     611   1831948030 :                t175 = t78*t83
     612   1831948030 :                t176 = t90**2
     613   1831948030 :                t177 = 0.1e1_dp/t176
     614   1831948030 :                t178 = t85*t177
     615   1831948030 :                t179 = A*t88
     616   1831948030 :                t182 = t83*t
     617   1831948030 :                t183 = t87*t182
     618   1831948030 :                t186 = t167 + t170 + 0.2e1_dp*t179*Arho + 0.4e1_dp*t183*trho
     619   1831948030 :                t187 = t178*t186
     620   1831948030 :                t189 = 0.2e1_dp*t162*t164 + t78*t83*t171*t91 - t175*t187
     621   1831948030 :                t190 = gamma_var*t189
     622   1831948030 :                t191 = 0.1e1_dp/t94
     623   1831948030 :                epsilon_cGGArho = e_c_u_0rho + t190*t191
     624   1831948030 :                kfrho = k_frho
     625   1831948030 :                ex_unifrho = -0.3e1_dp/0.4e1_dp*t5*kfrho
     626   1831948030 :                t195 = kf**2
     627   1831948030 :                t196 = 0.1e1_dp/t195
     628   1831948030 :                t197 = my_norm_drho*t196
     629   1831948030 :                t198 = t6*kfrho
     630   1831948030 :                t200 = t99*t115
     631   1831948030 :                srho = -t197*t198/0.2e1_dp - t200/0.2e1_dp
     632   1831948030 :                t202 = t105**2
     633   1831948030 :                t204 = 0.1e1_dp/t202*mu
     634   1831948030 :                Fxrho = 0.2e1_dp*t204*s*srho
     635   1831948030 :                t208 = my_rho*ex_unifrho
     636              : 
     637   1831948030 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     638              :                   e_rho(ii) = e_rho(ii) + &
     639              :                               scale_ex*(ex_unif*Fx + t208*Fx + t108*Fxrho) + &
     640   1701314678 :                               scale_ec*(epsilon_cGGA + my_rho*epsilon_cGGArho)
     641              :                END IF
     642              : 
     643   1831948030 :                tnorm_drho = t74*t6/0.2e1_dp
     644   1831948030 :                t214 = t163*tnorm_drho
     645   1831948030 :                t217 = t78*t182
     646   1831948030 :                t218 = A*tnorm_drho
     647   1831948030 :                t226 = 0.2e1_dp*t168*tnorm_drho + 0.4e1_dp*t183*tnorm_drho
     648              :                t229 = 0.2e1_dp*t162*t214 + 0.2e1_dp*t217*t218*t91 - &
     649   1831948030 :                       t175*t178*t226
     650   1831948030 :                t230 = gamma_var*t229
     651   1831948030 :                Hnorm_drho = t230*t191
     652   1831948030 :                snorm_drho = t98*t6/0.2e1_dp
     653   1831948030 :                Fxnorm_drho = 0.2e1_dp*t204*s*snorm_drho
     654              : 
     655   1831948030 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     656              :                   e_ndrho(ii) = e_ndrho(ii) + &
     657              :                                 scale_ex*t108*Fxnorm_drho + scale_ec*my_rho* &
     658   1701314678 :                                 Hnorm_drho
     659              :                END IF
     660              : 
     661   1831948030 :                t238 = 0.1e1_dp/t69
     662   1831948030 :                t239 = 0.1e1_dp/t111/t7*t238
     663   1831948030 :                t240 = t114**2
     664   1831948030 :                t241 = 0.1e1_dp/t240
     665   1831948030 :                t246 = 0.1e1_dp/t114/my_rho
     666              :                rsrhorho = -t4*t239*t241/0.18e2_dp + t4*t113* &
     667   1831948030 :                           t246/0.6e1_dp
     668   1831948030 :                t252 = 0.2e1_dp*t119*rsrhorho*t27
     669   1831948030 :                t253 = alpha_1_1*rsrho
     670   1831948030 :                t255 = t124*t138*t139
     671   1831948030 :                t259 = 0.1e1_dp/t123/t22
     672   1831948030 :                t260 = t11*t259
     673   1831948030 :                t261 = t138**2
     674   1831948030 :                t262 = t261*t139
     675   1831948030 :                t265 = 0.1e1_dp/t17
     676   1831948030 :                t266 = beta_1_1*t265
     677   1831948030 :                t267 = rsrho**2
     678   1831948030 :                t271 = t127*rsrhorho/0.2e1_dp
     679   1831948030 :                t272 = beta_2_1*rsrhorho
     680   1831948030 :                t273 = beta_3_1*t126
     681   1831948030 :                t277 = 0.3e1_dp/0.2e1_dp*t131*rsrhorho
     682   1831948030 :                t278 = t19**2
     683   1831948030 :                t280 = rs**2
     684   1831948030 :                t281 = 0.1e1_dp/t280
     685   1831948030 :                t286 = t21*t19*rsrhorho*t135
     686              :                t290 = -t266*t267/0.4e1_dp + t271 + t272 + 0.3e1_dp/0.4e1_dp &
     687              :                       *t273*t267 + t277 + t21*t278*t267*t281 + t286 - t21*t19 &
     688   1831948030 :                       *t267*t281
     689   1831948030 :                t291 = t290*t139
     690   1831948030 :                t293 = t123**2
     691   1831948030 :                t294 = 0.1e1_dp/t293
     692   1831948030 :                t295 = t11*t294
     693   1831948030 :                t296 = t26**2
     694   1831948030 :                t297 = 0.1e1_dp/t296
     695   1831948030 :                t299 = t261*t297*t13
     696              :                e_c_u_0rhorho = -t252 + 0.2e1_dp*t253*t255 - 0.2e1_dp*t260* &
     697   1831948030 :                                t262 + t125*t291 + t295*t299/0.2e1_dp
     698   1831948030 :                e_c_u_01rho = e_c_u_0rho
     699   1831948030 :                t305 = t69**2
     700   1831948030 :                k_frhorho = -0.2e1_dp/0.9e1_dp*t1/t142/t70*t305
     701   1831948030 :                t309 = 0.1e1_dp/t73/t72
     702   1831948030 :                t310 = k_frho**2
     703   1831948030 :                t315 = t146*k_frhorho*t5
     704   1831948030 :                k_srhorho = -t309*t310*t238/0.2e1_dp + t315
     705   1831948030 :                k_s1rho = k_srho
     706   1831948030 :                t317 = 0.1e1_dp/t148/k_s
     707   1831948030 :                t318 = my_norm_drho*t317
     708   1831948030 :                t321 = t115*k_srho
     709   1831948030 :                t323 = t150*t321/0.2e1_dp
     710   1831948030 :                t324 = t6*k_srhorho
     711   1831948030 :                t327 = t115*k_s1rho
     712   1831948030 :                t329 = t150*t327/0.2e1_dp
     713   1831948030 :                t330 = t75*t246
     714              :                trhorho = t318*t151*k_s1rho + t323 - t150*t324/0.2e1_dp + &
     715   1831948030 :                          t329 + t330
     716   1831948030 :                t331 = t6*k_s1rho
     717   1831948030 :                t1rho = -t150*t331/0.2e1_dp - t153/0.2e1_dp
     718   1831948030 :                t336 = beta/t155/gamma_var
     719   1831948030 :                t338 = 0.1e1_dp/t158/t81
     720   1831948030 :                t339 = t336*t338
     721   1831948030 :                t340 = t80**2
     722   1831948030 :                t341 = e_c_u_0rho*t340
     723   1831948030 :                t348 = t336*t159
     724   1831948030 :                t349 = e_c_u_0rho*e_c_u_01rho
     725              :                Arhorho = 0.2e1_dp*t339*t341*e_c_u_01rho + t157*t159* &
     726   1831948030 :                          e_c_u_0rhorho*t80 - t348*t349*t80
     727   1831948030 :                A1rho = t157*t159*e_c_u_01rho*t80
     728   1831948030 :                t354 = t78*t1rho
     729   1831948030 :                t357 = A1rho*t83
     730   1831948030 :                t359 = 0.2e1_dp*t168*t1rho
     731   1831948030 :                t360 = t357 + t359
     732   1831948030 :                t361 = t360*t91
     733   1831948030 :                t362 = t361*trho
     734   1831948030 :                t369 = t357 + t359 + 0.2e1_dp*t179*A1rho + 0.4e1_dp*t183*t1rho
     735   1831948030 :                t370 = trho*t369
     736   1831948030 :                t371 = t178*t370
     737   1831948030 :                t374 = t163*trhorho
     738   1831948030 :                t377 = t171*t91
     739   1831948030 :                t378 = t377*t1rho
     740   1831948030 :                t381 = Arhorho*t83
     741   1831948030 :                t382 = Arho*t
     742   1831948030 :                t384 = 0.2e1_dp*t382*t1rho
     743   1831948030 :                t385 = A1rho*t
     744   1831948030 :                t387 = 0.2e1_dp*t385*trho
     745   1831948030 :                t388 = A*t1rho
     746   1831948030 :                t390 = 0.2e1_dp*t388*trho
     747   1831948030 :                t392 = 0.2e1_dp*t168*trhorho
     748   1831948030 :                t393 = t381 + t384 + t387 + t390 + t392
     749   1831948030 :                t397 = t171*t177
     750   1831948030 :                t400 = t186*t1rho
     751   1831948030 :                t401 = t178*t400
     752   1831948030 :                t404 = t360*t177
     753   1831948030 :                t408 = 0.1e1_dp/t176/t90
     754   1831948030 :                t409 = t85*t408
     755   1831948030 :                t410 = t186*t369
     756   1831948030 :                t414 = A1rho*t88
     757   1831948030 :                t417 = A*t182
     758   1831948030 :                t418 = Arho*t1rho
     759   1831948030 :                t423 = trho*A1rho
     760   1831948030 :                t426 = t87*t83
     761   1831948030 :                t427 = trho*t1rho
     762              :                t432 = t381 + t384 + t387 + t390 + t392 + 0.2e1_dp*t414*Arho + &
     763              :                       0.8e1_dp*t417*t418 + 0.2e1_dp*t179*Arhorho + 0.8e1_dp* &
     764   1831948030 :                       t417*t423 + 0.12e2_dp*t426*t427 + 0.4e1_dp*t183*trhorho
     765              :                t435 = 0.2e1_dp*t354*t164 + 0.2e1_dp*t162*t362 - 0.2e1_dp &
     766              :                       *t162*t371 + 0.2e1_dp*t162*t374 + 0.2e1_dp*t162*t378 + &
     767              :                       t78*t83*t393*t91 - t175*t397*t369 - 0.2e1_dp*t162*t401 &
     768              :                       - t175*t404*t186 + 0.2e1_dp*t175*t409*t410 - t175*t178 &
     769   1831948030 :                       *t432
     770   1831948030 :                t436 = gamma_var*t435
     771   1831948030 :                t438 = t94**2
     772   1831948030 :                t439 = 0.1e1_dp/t438
     773   1831948030 :                t440 = t163*t1rho
     774              :                t448 = 0.2e1_dp*t162*t440 + t78*t83*t360*t91 - t175* &
     775   1831948030 :                       t178*t369
     776   1831948030 :                t449 = t439*t448
     777   1831948030 :                t451 = gamma_var*t448
     778   1831948030 :                epsilon_cGGArhorho = e_c_u_0rhorho + t436*t191 - t190*t449
     779   1831948030 :                kfrhorho = k_frhorho
     780   1831948030 :                ex_unifrhorho = -0.3e1_dp/0.4e1_dp*t5*kfrhorho
     781   1831948030 :                ex_unif1rho = ex_unifrho
     782   1831948030 :                t456 = 0.1e1_dp/t195/kf
     783   1831948030 :                t457 = my_norm_drho*t456
     784   1831948030 :                t458 = kfrho**2
     785   1831948030 :                t459 = t6*t458
     786   1831948030 :                t461 = t115*kfrho
     787   1831948030 :                t462 = t197*t461
     788   1831948030 :                t463 = t6*kfrhorho
     789   1831948030 :                t465 = t197*t463/0.2e1_dp
     790   1831948030 :                t466 = t99*t246
     791   1831948030 :                srhorho = t457*t459 + t462 - t465 + t466
     792   1831948030 :                s1rho = srho
     793   1831948030 :                t469 = mu**2
     794   1831948030 :                t470 = 0.1e1_dp/t202/t105*t469
     795   1831948030 :                t471 = t470*t101
     796   1831948030 :                t472 = srho*t103
     797   1831948030 :                t476 = s1rho*srho
     798              :                Fxrhorho = -0.8e1_dp*t471*t472*s1rho + 0.2e1_dp*t204* &
     799   1831948030 :                           t476 + 0.2e1_dp*t204*s*srhorho
     800   1831948030 :                Fx1rho = 0.2e1_dp*t204*s*s1rho
     801   1831948030 :                t487 = my_rho*ex_unifrhorho
     802   1831948030 :                t491 = my_rho*ex_unif1rho
     803              : 
     804   1831948030 :                IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     805              :                   e_rho_rho(ii) = e_rho_rho(ii) + &
     806              :                                   scale_ex*(ex_unif1rho*Fx + ex_unif*Fx1rho + &
     807              :                                             ex_unifrho*Fx + t487*Fx + t208*Fx1rho + ex_unif*Fxrho + t491 &
     808              :                                             *Fxrho + t108*Fxrhorho) + scale_ec*(e_c_u_01rho + t451*t191 &
     809    187411902 :                                                                                 + epsilon_cGGArho + my_rho*epsilon_cGGArhorho)
     810              :                END IF
     811              : 
     812   1831948030 :                t496 = t149*t6
     813   1831948030 :                t498 = t74*t115
     814   1831948030 :                tnorm_drhorho = -t496*k_srho/0.2e1_dp - t498/0.2e1_dp
     815   1831948030 :                t500 = t78*trho
     816   1831948030 :                t503 = t377*tnorm_drho
     817   1831948030 :                t506 = tnorm_drho*t186
     818   1831948030 :                t507 = t178*t506
     819   1831948030 :                t510 = t163*tnorm_drhorho
     820   1831948030 :                t513 = t91*trho
     821   1831948030 :                t517 = Arho*tnorm_drho
     822   1831948030 :                t521 = A*tnorm_drhorho
     823   1831948030 :                t525 = t177*t186
     824   1831948030 :                t529 = t226*trho
     825   1831948030 :                t530 = t178*t529
     826   1831948030 :                t535 = t226*t186
     827   1831948030 :                t541 = A*trho
     828   1831948030 :                t548 = tnorm_drho*trho
     829              :                t553 = 0.2e1_dp*t382*tnorm_drho + 0.2e1_dp*t541*tnorm_drho &
     830              :                       + 0.2e1_dp*t168*tnorm_drhorho + 0.8e1_dp*t417*t517 + &
     831   1831948030 :                       0.12e2_dp*t426*t548 + 0.4e1_dp*t183*tnorm_drhorho
     832              :                t556 = 0.2e1_dp*t500*t214 + 0.2e1_dp*t162*t503 - 0.2e1_dp &
     833              :                       *t162*t507 + 0.2e1_dp*t162*t510 + 0.6e1_dp*t175*t218* &
     834              :                       t513 + 0.2e1_dp*t217*t517*t91 + 0.2e1_dp*t217*t521*t91 - &
     835              :                       0.2e1_dp*t217*t218*t525 - 0.2e1_dp*t162*t530 - t175* &
     836   1831948030 :                       t397*t226 + 0.2e1_dp*t175*t409*t535 - t175*t178*t553
     837   1831948030 :                t557 = gamma_var*t556
     838   1831948030 :                t559 = t439*t189
     839   1831948030 :                Hnorm_drhorho = t557*t191 - t230*t559
     840   1831948030 :                t562 = t196*t6
     841   1831948030 :                snorm_drhorho = -t562*kfrho/0.2e1_dp - t98*t115/0.2e1_dp
     842   1831948030 :                t566 = snorm_drho*t103
     843              :                Fxnorm_drhorho = -0.8e1_dp*t471*t566*srho + 0.2e1_dp*t204 &
     844   1831948030 :                                 *srho*snorm_drho + 0.2e1_dp*t204*s*snorm_drhorho
     845              : 
     846   1831948030 :                IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     847              :                   e_ndrho_rho(ii) = e_ndrho_rho(ii) + &
     848              :                                     scale_ex*(ex_unif*Fxnorm_drho + t208* &
     849              :                                               Fxnorm_drho + t108*Fxnorm_drhorho) + scale_ec*(Hnorm_drho + my_rho &
     850    187411902 :                                                                                              *Hnorm_drhorho)
     851              :                END IF
     852              : 
     853   1831948030 :                t581 = tnorm_drho**2
     854   1831948030 :                t586 = A*t581
     855   1831948030 :                t590 = tnorm_drho*t226
     856   1831948030 :                t591 = t178*t590
     857   1831948030 :                t594 = t177*t226
     858   1831948030 :                t598 = t226**2
     859   1831948030 :                t605 = 0.2e1_dp*t586 + 0.12e2_dp*t426*t581
     860              :                t609 = gamma_var*(0.2e1_dp*t78*t581*t85*t91 + 0.10e2_dp &
     861              :                                  *t175*t586*t91 - 0.4e1_dp*t162*t591 - 0.4e1_dp*t217* &
     862   1831948030 :                                  t218*t594 + 0.2e1_dp*t175*t409*t598 - t175*t178*t605)
     863   1831948030 :                t611 = t229**2
     864   1831948030 :                t612 = gamma_var*t611
     865   1831948030 :                Hnorm_drhonorm_drho = t609*t191 - t612*t439
     866   1831948030 :                t614 = snorm_drho**2
     867              :                Fxnorm_drhonorm_drho = -0.8e1_dp*t470*t101*t614*t103 + &
     868   1831948030 :                                       0.2e1_dp*t204*t614
     869              : 
     870   1831948030 :                IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     871              :                   e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
     872              :                                       scale_ex*t108*Fxnorm_drhonorm_drho + &
     873    187411902 :                                       scale_ec*my_rho*Hnorm_drhonorm_drho
     874              :                END IF
     875              : 
     876   1831948030 :                IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
     877              :                   rsrhorhorho = -0.5e1_dp/0.54e2_dp*t4/t111/t238/ &
     878              :                                 t115*t627/t240/t114 + t4*t239/t240/my_rho/0.3e1_dp &
     879            0 :                                 - t4*t113*t241/0.2e1_dp
     880            0 :                   rs2rho = rsrho
     881            0 :                   t645 = alpha_1_1*rsrhorho
     882              :                   t654 = t127*rs2rho/0.2e1_dp + beta_2_1*rs2rho + 0.3e1_dp/ &
     883            0 :                          0.2e1_dp*t131*rs2rho + t21*t19*rs2rho*t135
     884            0 :                   t656 = t124*t654*t139
     885            0 :                   t661 = t140*t654
     886            0 :                   t664 = rsrho*rs2rho
     887            0 :                   t669 = t21*t278
     888            0 :                   t670 = rs2rho*t281
     889            0 :                   t671 = t670*rsrho
     890            0 :                   t673 = t21*t19
     891              :                   t675 = -t266*t664/0.4e1_dp + t271 + t272 + 0.3e1_dp/0.4e1_dp &
     892            0 :                          *t273*t664 + t277 + t669*t671 + t286 - t673*t671
     893            0 :                   t685 = alpha_1_1*rs2rho
     894            0 :                   t693 = t675*t139
     895            0 :                   t701 = t297*t13
     896            0 :                   t702 = t701*t654
     897            0 :                   t714 = t267*rs2rho
     898            0 :                   t717 = rsrhorho*rsrho
     899            0 :                   t720 = rsrhorho*rs2rho
     900            0 :                   t740 = rs2rho/t280/rs*t267
     901            0 :                   t743 = rsrhorho*t281*rsrho
     902            0 :                   t748 = t670*rsrhorho
     903              :                   t758 = 0.3e1_dp/0.8e1_dp*beta_1_1/t14/t280*t714 - t266* &
     904              :                          t717/0.2e1_dp - t266*t720/0.4e1_dp + t127*rsrhorhorho/ &
     905              :                          0.2e1_dp + beta_2_1*rsrhorhorho - 0.3e1_dp/0.8e1_dp*beta_3_1* &
     906              :                          t265*t714 + 0.3e1_dp/0.2e1_dp*t273*t717 + 0.3e1_dp/ &
     907              :                          0.4e1_dp*t273*t720 + 0.3e1_dp/0.2e1_dp*t131*rsrhorhorho + &
     908              :                          t21*t278*t19*t740 + 0.2e1_dp*t669*t743 - 0.3e1_dp*t669* &
     909              :                          t740 + t669*t748 + t21*t19*rsrhorhorho*t135 - t673*t748 - &
     910            0 :                          0.2e1_dp*t673*t743 + 0.2e1_dp*t673*t740
     911            0 :                   t776 = A_1**2
     912              :                   e_c_u_0rhorhorho = -0.2e1_dp*t119*rsrhorhorho*t27 + t645* &
     913              :                                      t656 + 0.2e1_dp*t645*t255 - 0.4e1_dp*t253*t259*t661 + &
     914              :                                      0.2e1_dp*t253*t124*t675*t139 + t253*t294*t138*t297* &
     915              :                                      t13*t654 - 0.2e1_dp*t685*t259*t261*t139 + 0.6e1_dp*t295 &
     916              :                                      *t262*t654 - 0.4e1_dp*t260*t693*t138 - 0.3e1_dp*t11/ &
     917              :                                      t293/t22*t261*t702 + t685*t124*t290*t139 - 0.2e1_dp* &
     918              :                                      t260*t291*t654 + t125*t758*t139 + t295*t290*t702/ &
     919              :                                      0.2e1_dp + t685*t294*t299/0.2e1_dp + t295*t675*t701*t138 &
     920            0 :                                      + t11/t293/t123*t261/t296/t26/t776*t654/0.2e1_dp
     921              :                   e_c_u_0rho1rho = -t252 + t253*t656 + t685*t255 - 0.2e1_dp* &
     922            0 :                                    t260*t661 + t125*t693 + t295*t138*t702/0.2e1_dp
     923            0 :                   e_c_u_01rhorho = e_c_u_0rho1rho
     924            0 :                   e_c_u_02rho = -0.2e1_dp*t119*rs2rho*t27 + t125*t654*t139
     925            0 :                   k_frhorhorho = 0.10e2_dp/0.27e2_dp*t1/t142/t114*t69
     926            0 :                   k_f2rho = kfrho
     927            0 :                   t801 = k_f**2
     928            0 :                   t809 = t309*k_frhorho
     929            0 :                   t812 = t238*k_f2rho
     930            0 :                   k_srho1rho = -t309*k_frho*t812/0.2e1_dp + t315
     931            0 :                   k_s1rhorho = k_srho1rho
     932            0 :                   k_s2rho = t146*k_f2rho*t5
     933            0 :                   t821 = t148**2
     934            0 :                   t825 = k_srho*k_s1rho
     935            0 :                   t831 = t6*k_srho1rho
     936              :                   trhorhorho = -0.3e1_dp*my_norm_drho/t821*t6*t825*k_s2rho - &
     937              :                                t318*t321*k_s1rho + t318*t831*k_s1rho + t318*t151* &
     938              :                                k_s1rhorho - t318*t321*k_s2rho - t150*t246*k_srho + t150* &
     939              :                                t115*k_srho1rho/0.2e1_dp + t318*t324*k_s2rho + t150*t115* &
     940              :                                k_srhorho/0.2e1_dp - t150*t6*(0.3e1_dp/0.4e1_dp/t73/ &
     941              :                                                              t801/t238*t310*t627*k_f2rho - t809*t238*k_frho - t809* &
     942              :                                                              t812/0.2e1_dp + t146*k_frhorhorho*t5)/0.2e1_dp - t318*t327 &
     943              :                                *k_s2rho - t150*t246*k_s1rho + t150*t115*k_s1rhorho/ &
     944            0 :                                0.2e1_dp - t150*t246*k_s2rho - 0.3e1_dp*t75*t241
     945            0 :                   t868 = t150*t115*k_s2rho/0.2e1_dp
     946              :                   trho1rho = t318*t151*k_s2rho + t323 - t150*t831/0.2e1_dp + &
     947            0 :                              t868 + t330
     948              :                   t1rhorho = t318*t331*k_s2rho + t329 - t150*t6*k_s1rhorho/ &
     949            0 :                              0.2e1_dp + t868 + t330
     950            0 :                   t2rho = -t150*t6*k_s2rho/0.2e1_dp - t153/0.2e1_dp
     951            0 :                   t877 = t155**2
     952            0 :                   t879 = beta/t877
     953            0 :                   t880 = t158**2
     954            0 :                   t885 = e_c_u_01rho*e_c_u_02rho
     955              :                   Arhorhorho = 0.6e1_dp*t879/t880*e_c_u_0rho*t340*t80* &
     956              :                                t885 + 0.2e1_dp*t339*e_c_u_0rho1rho*t340*e_c_u_01rho - &
     957              :                                0.6e1_dp*t879*t338*t341*t885 + 0.2e1_dp*t339*t341* &
     958              :                                e_c_u_01rhorho + 0.2e1_dp*t339*e_c_u_0rhorho*t340* &
     959              :                                e_c_u_02rho + t157*t159*e_c_u_0rhorhorho*t80 - t348* &
     960              :                                e_c_u_0rhorho*e_c_u_02rho*t80 - t348*e_c_u_0rho1rho* &
     961              :                                e_c_u_01rho*t80 - t348*e_c_u_0rho*e_c_u_01rhorho*t80 + t879 &
     962            0 :                                *t159*t349*e_c_u_02rho*t80
     963              :                   Arho1rho = 0.2e1_dp*t339*t341*e_c_u_02rho + t157*t159* &
     964            0 :                              e_c_u_0rho1rho*t80 - t348*e_c_u_0rho*e_c_u_02rho*t80
     965              :                   A1rhorho = 0.2e1_dp*t339*e_c_u_01rho*t340*e_c_u_02rho + &
     966            0 :                              t157*t159*e_c_u_01rhorho*t80 - t348*t885*t80
     967            0 :                   A2rho = t157*t159*e_c_u_02rho*t80
     968            0 :                   t940 = Arho1rho*t83
     969            0 :                   t942 = 0.2e1_dp*t382*t2rho
     970            0 :                   t943 = A2rho*t
     971            0 :                   t945 = 0.2e1_dp*t943*trho
     972            0 :                   t946 = A*t2rho
     973            0 :                   t948 = 0.2e1_dp*t946*trho
     974            0 :                   t950 = 0.2e1_dp*t168*trho1rho
     975            0 :                   t951 = A2rho*t88
     976            0 :                   t954 = Arho*t2rho
     977              :                   t967 = t940 + t942 + t945 + t948 + t950 + 0.2e1_dp*t951*Arho + &
     978              :                          0.8e1_dp*t417*t954 + 0.2e1_dp*t179*Arho1rho + 0.8e1_dp* &
     979              :                          t417*trho*A2rho + 0.12e2_dp*t426*trho*t2rho + 0.4e1_dp* &
     980            0 :                          t183*trho1rho
     981            0 :                   t976 = t78*t2rho
     982            0 :                   t980 = t78*t*t85
     983            0 :                   t982 = A2rho*t83
     984            0 :                   t984 = 0.2e1_dp*t168*t2rho
     985            0 :                   t989 = t982 + t984 + 0.2e1_dp*t179*A2rho + 0.4e1_dp*t183*t2rho
     986            0 :                   t990 = t369*t989
     987            0 :                   t994 = t982 + t984
     988            0 :                   t995 = t994*t177
     989            0 :                   t998 = Arhorhorho*t83
     990            0 :                   t999 = Arhorho*t
     991            0 :                   t1001 = 0.2e1_dp*t999*t2rho
     992            0 :                   t1004 = 0.2e1_dp*Arho1rho*t*t1rho
     993            0 :                   t1005 = t954*t1rho
     994            0 :                   t1006 = 0.2e1_dp*t1005
     995            0 :                   t1008 = 0.2e1_dp*t382*t1rhorho
     996            0 :                   t1011 = 0.2e1_dp*A1rhorho*t*trho
     997            0 :                   t1012 = A1rho*t2rho
     998            0 :                   t1013 = t1012*trho
     999            0 :                   t1014 = 0.2e1_dp*t1013
    1000            0 :                   t1016 = 0.2e1_dp*t385*trho1rho
    1001            0 :                   t1017 = A2rho*t1rho
    1002            0 :                   t1018 = t1017*trho
    1003            0 :                   t1019 = 0.2e1_dp*t1018
    1004            0 :                   t1022 = 0.2e1_dp*A*t1rhorho*trho
    1005            0 :                   t1024 = 0.2e1_dp*t388*trho1rho
    1006            0 :                   t1026 = 0.2e1_dp*t943*trhorho
    1007            0 :                   t1028 = 0.2e1_dp*t946*trhorho
    1008            0 :                   t1030 = 0.2e1_dp*t168*trhorhorho
    1009              :                   t1031 = t998 + t1001 + t1004 + t1006 + t1008 + t1011 + t1014 + &
    1010            0 :                           t1016 + t1019 + t1022 + t1024 + t1026 + t1028 + t1030
    1011            0 :                   t1035 = t940 + t942 + t945 + t948 + t950
    1012            0 :                   t1042 = t369*t2rho
    1013            0 :                   t1046 = A1rhorho*t83
    1014            0 :                   t1048 = 0.2e1_dp*t385*t2rho
    1015            0 :                   t1050 = 0.2e1_dp*t943*t1rho
    1016            0 :                   t1052 = 0.2e1_dp*t946*t1rho
    1017            0 :                   t1054 = 0.2e1_dp*t168*t1rhorho
    1018            0 :                   t1055 = t1046 + t1048 + t1050 + t1052 + t1054
    1019            0 :                   t1060 = t78*t86
    1020            0 :                   t1061 = t176**2
    1021            0 :                   t1062 = 0.1e1_dp/t1061
    1022              :                   t1067 = -0.2e1_dp*t162*t178*t967*t1rho + 0.2e1_dp*t175* &
    1023              :                           t409*t967*t369 + 0.2e1_dp*t976*t374 + 0.4e1_dp*t980* &
    1024              :                           t408*trho*t990 - t175*t995*t432 + t78*t83*t1031*t91 - &
    1025              :                           t175*t1035*t177*t369 - 0.2e1_dp*t162*t995*t370 - &
    1026              :                           0.2e1_dp*t162*t397*t1042 + 0.2e1_dp*t162*t1055*t91* &
    1027            0 :                           trho - 0.6e1_dp*t1060*t1062*t186*t990
    1028              :                   t1093 = t1046 + t1048 + t1050 + t1052 + t1054 + 0.2e1_dp*t951* &
    1029              :                           A1rho + 0.8e1_dp*t417*t1012 + 0.2e1_dp*t179*A1rhorho + &
    1030              :                           0.8e1_dp*t417*t1017 + 0.12e2_dp*t426*t1rho*t2rho + &
    1031            0 :                           0.4e1_dp*t183*t1rhorho
    1032            0 :                   t1097 = t1rho*t989
    1033            0 :                   t1103 = trho*t989
    1034            0 :                   t1104 = t178*t1103
    1035            0 :                   t1106 = t994*t91
    1036            0 :                   t1109 = t171*t408
    1037              :                   t1115 = t976*t362 + t162*t361*trho1rho - t162*t178*t432 &
    1038              :                           *t2rho + t175*t994*t408*t410 - t162*t178*trho1rho*t369 &
    1039              :                           - t162*t178*trho*t1093 - t162*t397*t1097 + t175*t409* &
    1040              :                           t186*t1093 - t354*t1104 + t162*t1106*trhorho + t175*t1109 &
    1041            0 :                           *t990 + t175*t409*t432*t989
    1042            0 :                   t1118 = t1106*trho
    1043            0 :                   t1121 = t360*t408
    1044            0 :                   t1122 = t186*t989
    1045            0 :                   t1126 = t393*t177
    1046            0 :                   t1129 = t408*t186
    1047            0 :                   t1148 = t186*t2rho
    1048              :                   t1152 = 0.2e1_dp*t354*t1118 + 0.2e1_dp*t175*t1121*t1122 &
    1049              :                           - t175*t1126*t989 + 0.4e1_dp*t980*t1129*t1042 + 0.2e1_dp* &
    1050              :                           t976*t378 - t175*t404*t967 + 0.2e1_dp*t162*t377* &
    1051              :                           t1rhorho - 0.2e1_dp*t976*t401 + 0.2e1_dp*t78*t1rhorho*t164 &
    1052            0 :                           - 0.2e1_dp*t162*t404*t1103 - 0.2e1_dp*t162*t404*t1148
    1053            0 :                   t1157 = t393*t91
    1054            0 :                   t1167 = A2rho*t182
    1055            0 :                   t1187 = A1rho*t182
    1056            0 :                   t1196 = t87*t
    1057              :                   t1203 = t1008 + 0.8e1_dp*t417*trhorho*A2rho + 0.12e2_dp* &
    1058              :                           t426*trhorho*t2rho + 0.8e1_dp*t1167*t423 + 0.8e1_dp*t417* &
    1059              :                           trho*A1rhorho + 0.8e1_dp*t417*Arho*t1rhorho + 0.8e1_dp* &
    1060              :                           t1167*t418 + 0.8e1_dp*t417*Arho1rho*t1rho + 0.12e2_dp*t426 &
    1061              :                           *trho1rho*t1rho + 0.12e2_dp*t426*trho*t1rhorho + t998 + &
    1062              :                           0.8e1_dp*t1187*t954 + 0.8e1_dp*t417*trho1rho*A1rho + &
    1063              :                           0.8e1_dp*t417*Arhorho*t2rho + 0.24e2_dp*t1196*t427*t2rho &
    1064            0 :                           + 0.2e1_dp*A1rhorho*t88*Arho + t1014
    1065              :                   t1218 = t1016 + t1001 + 0.2e1_dp*t951*Arhorho + t1026 + t1028 &
    1066              :                           + t1022 + t1011 + t1030 + 0.2e1_dp*t179*Arhorhorho + 0.4e1_dp* &
    1067              :                           t183*trhorhorho + 0.2e1_dp*t414*Arho1rho + t1004 + t1006 + &
    1068              :                           t1019 + t1024 + 0.24e2_dp*t84*t1018 + 0.24e2_dp*t84*t1005 + &
    1069            0 :                           0.24e2_dp*t84*t1013
    1070            0 :                   t1226 = t163*trho1rho
    1071              :                   t1249 = -0.2e1_dp*t162*t178*t186*t1rhorho + 0.2e1_dp* &
    1072              :                           t162*t1157*t2rho - t175*t178*(t1203 + t1218) + 0.2e1_dp* &
    1073              :                           t162*t1035*t91*t1rho + 0.2e1_dp*t354*t1226 - 0.2e1_dp* &
    1074              :                           t162*t995*t400 + 0.2e1_dp*t162*t163*trhorhorho - 0.2e1_dp &
    1075              :                           *t162*t178*trhorho*t989 - t175*t1055*t177*t186 - &
    1076              :                           0.2e1_dp*t976*t371 - t175*t397*t1093 + 0.4e1_dp*t980* &
    1077            0 :                           t1129*t1097
    1078              :                   t1262 = 0.2e1_dp*t162*t163*t2rho + t78*t83*t994*t91 - &
    1079            0 :                           t175*t178*t989
    1080            0 :                   t1263 = t439*t1262
    1081              :                   t1291 = 0.2e1_dp*t976*t164 + 0.2e1_dp*t162*t1118 - &
    1082              :                           0.2e1_dp*t162*t1104 + 0.2e1_dp*t162*t1226 + 0.2e1_dp*t162 &
    1083              :                           *t377*t2rho + t78*t83*t1035*t91 - t175*t397*t989 - &
    1084              :                           0.2e1_dp*t162*t178*t1148 - t175*t995*t186 + 0.2e1_dp* &
    1085            0 :                           t175*t409*t1122 - t175*t178*t967
    1086            0 :                   t1292 = gamma_var*t1291
    1087            0 :                   t1295 = 0.1e1_dp/t438/t94
    1088              :                   t1329 = 0.2e1_dp*t976*t440 + 0.2e1_dp*t162*t1106*t1rho - &
    1089              :                           0.2e1_dp*t162*t178*t1097 + 0.2e1_dp*t162*t163*t1rhorho &
    1090              :                           + 0.2e1_dp*t162*t361*t2rho + t78*t83*t1055*t91 - t175* &
    1091              :                           t404*t989 - 0.2e1_dp*t162*t178*t1042 - t175*t995*t369 + &
    1092            0 :                           0.2e1_dp*t175*t409*t990 - t175*t178*t1093
    1093            0 :                   kfrhorhorho = k_frhorhorho
    1094            0 :                   kf2rho = k_f2rho
    1095            0 :                   ex_unifrho1rho = ex_unifrhorho
    1096            0 :                   ex_unif1rhorho = ex_unifrho1rho
    1097            0 :                   ex_unif2rho = -0.3e1_dp/0.4e1_dp*t5*kf2rho
    1098            0 :                   t1342 = t195**2
    1099              :                   srho1rho = t457*t198*kf2rho + t462/0.2e1_dp - t465 + t197* &
    1100            0 :                              t115*kf2rho/0.2e1_dp + t466
    1101            0 :                   s1rhorho = srho1rho
    1102            0 :                   s2rho = -t197*t6*kf2rho/0.2e1_dp - t200/0.2e1_dp
    1103            0 :                   t1380 = t202**2
    1104            0 :                   t1385 = 0.1e1_dp/t1380*t469*mu*t101*s
    1105            0 :                   t1386 = kappa**2
    1106            0 :                   t1387 = 0.1e1_dp/t1386
    1107            0 :                   t1389 = s1rho*s2rho
    1108            0 :                   t1393 = t470*s
    1109              :                   Fxrho1rho = -0.8e1_dp*t471*t472*s2rho + 0.2e1_dp*t204* &
    1110            0 :                               s2rho*srho + 0.2e1_dp*t204*s*srho1rho
    1111              :                   Fx1rhorho = -0.8e1_dp*t471*s1rho*t103*s2rho + 0.2e1_dp* &
    1112            0 :                               t204*t1389 + 0.2e1_dp*t204*s*s1rhorho
    1113            0 :                   Fx2rho = 0.2e1_dp*t204*s*s2rho
    1114              :                   ex_ldarhorhorho = ex_unif1rhorho*Fx + ex_unif1rho*Fx2rho + &
    1115              :                                     ex_unif2rho*Fx1rho + ex_unif*Fx1rhorho + ex_unifrho1rho*Fx + &
    1116              :                                     ex_unifrho*Fx2rho + ex_unifrhorho*Fx - 0.3e1_dp/0.4e1_dp*my_rho &
    1117              :                                     *t5*kfrhorhorho*Fx + t487*Fx2rho + ex_unifrho*Fx1rho + my_rho &
    1118              :                                     *ex_unifrho1rho*Fx1rho + t208*Fx1rhorho + ex_unif2rho*Fxrho &
    1119              :                                     + ex_unif*Fxrho1rho + ex_unif1rho*Fxrho + my_rho*ex_unif1rhorho* &
    1120              :                                     Fxrho + t491*Fxrho1rho + ex_unif*Fxrhorho + my_rho*ex_unif2rho* &
    1121              :                                     Fxrhorho + t108*(0.48e2_dp*t1385*srho*t1387*t1389 - &
    1122              :                                                      0.24e2_dp*t1393*t472*t1389 - 0.8e1_dp*t471*srho1rho*t103 &
    1123              :                                                      *s1rho - 0.8e1_dp*t471*t472*s1rhorho + 0.2e1_dp*t204* &
    1124              :                                                      s1rhorho*srho + 0.2e1_dp*t204*s1rho*srho1rho - 0.8e1_dp* &
    1125              :                                                      t471*srhorho*t103*s2rho + 0.2e1_dp*t204*s2rho*srhorho + &
    1126              :                                                      0.2e1_dp*t204*s*(-0.3e1_dp*my_norm_drho/t1342*t459*kf2rho &
    1127              :                                                                       - t457*t115*t458 + 0.2e1_dp*t457*t463*kfrho - 0.2e1_dp* &
    1128              :                                                                       t457*t461*kf2rho - 0.2e1_dp*t197*t246*kfrho + 0.3e1_dp/ &
    1129              :                                                                       0.2e1_dp*t197*t115*kfrhorho + t457*t463*kf2rho - t197*t6 &
    1130              :                                                                       *kfrhorhorho/0.2e1_dp - t197*t246*kf2rho - 0.3e1_dp*t99* &
    1131            0 :                                                                       t241))
    1132              : 
    1133              :                   e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + &
    1134              :                                       scale_ex*ex_ldarhorhorho + scale_ec*( &
    1135              :                                       e_c_u_01rhorho + gamma_var*t1329*t191 - t451*t1263 + &
    1136              :                                       e_c_u_0rho1rho + t1292*t191 - t190*t1263 + epsilon_cGGArhorho + &
    1137              :                                       my_rho*(e_c_u_0rhorhorho + gamma_var*(t1067 + 0.2e1_dp*t1115 + &
    1138              :                                                                          t1152 + t1249)*t191 - t436*t1263 - t1292*t449 + 0.2e1_dp* &
    1139            0 :                                               t190*t1295*t448*t1262 - t190*t439*t1329))
    1140              : 
    1141            0 :                   t1468 = t149*t115
    1142              :                   tnorm_drhorhorho = t317*t6*t825 + t1468*k_srho/0.2e1_dp - &
    1143              :                                      t496*k_srhorho/0.2e1_dp + t1468*k_s1rho/0.2e1_dp + t74* &
    1144            0 :                                      t246
    1145            0 :                   tnorm_drho1rho = -t496*k_s1rho/0.2e1_dp - t498/0.2e1_dp
    1146            0 :                   t1482 = A1rho*tnorm_drhorho
    1147            0 :                   t1492 = t78*t84
    1148            0 :                   t1493 = tnorm_drho*t177
    1149            0 :                   t1504 = t408*tnorm_drho
    1150            0 :                   t1505 = t1504*t410
    1151            0 :                   t1511 = A1rho*tnorm_drho
    1152            0 :                   t1515 = t226*t1rho
    1153            0 :                   t1521 = A*tnorm_drho1rho
    1154              :                   t1525 = 0.2e1_dp*t175*t409*t553*t369 + 0.2e1_dp*t217* &
    1155              :                           t1482*t91 - 0.2e1_dp*t162*t178*tnorm_drhorho*t369 + &
    1156              :                           0.2e1_dp*t354*t510 - 0.6e1_dp*t1492*t1493*t400 + 0.6e1_dp &
    1157              :                           *t175*t218*t91*trhorho + 0.2e1_dp*t162*t1157*tnorm_drho &
    1158              :                           + 0.4e1_dp*t980*t1505 - 0.6e1_dp*t1492*t1493*t370 + &
    1159              :                           0.6e1_dp*t175*t1511*t513 - 0.2e1_dp*t162*t397*t1515 - &
    1160            0 :                           0.2e1_dp*t354*t507 + 0.6e1_dp*t175*t1521*t513
    1161            0 :                   t1528 = Arhorho*tnorm_drho
    1162            0 :                   t1532 = t177*t369
    1163            0 :                   t1545 = t78*t417
    1164              :                   t1565 = 0.2e1_dp*t385*tnorm_drho + 0.2e1_dp*t388* &
    1165              :                           tnorm_drho + 0.2e1_dp*t168*tnorm_drho1rho + 0.8e1_dp*t417* &
    1166              :                           t1511 + 0.12e2_dp*t426*tnorm_drho*t1rho + 0.4e1_dp*t183* &
    1167            0 :                           tnorm_drho1rho
    1168            0 :                   t1573 = t91*t1rho
    1169              :                   t1584 = -0.2e1_dp*t354*t530 + 0.2e1_dp*t217*t1528*t91 - &
    1170              :                           0.2e1_dp*t217*t517*t1532 - 0.2e1_dp*t217*t1511*t525 + &
    1171              :                           0.2e1_dp*t162*t377*tnorm_drho1rho + 0.2e1_dp*t162*t361* &
    1172              :                           tnorm_drhorho + 0.4e1_dp*t1545*t1505 + 0.2e1_dp*t217*A* &
    1173              :                           tnorm_drhorhorho*t91 + 0.2e1_dp*t175*t409*t1565*t186 - &
    1174              :                           0.2e1_dp*t217*t521*t1532 + 0.6e1_dp*t175*t521*t1573 - &
    1175              :                           0.6e1_dp*t1060*t1062*t226*t410 + 0.6e1_dp*t175*t517* &
    1176            0 :                           t1573
    1177            0 :                   t1608 = t408*t226
    1178            0 :                   t1612 = t163*tnorm_drho1rho
    1179              :                   t1628 = -0.2e1_dp*t162*t404*t506 + 0.2e1_dp*t354*t503 - &
    1180              :                           0.2e1_dp*t162*t178*tnorm_drho*t432 - 0.2e1_dp*t162*t178 &
    1181              :                           *t553*t1rho - 0.2e1_dp*t162*t404*t529 - 0.2e1_dp*t162* &
    1182              :                           t178*t1565*trho - t175*t1126*t226 + 0.4e1_dp*t980*t1608 &
    1183              :                           *t370 + 0.2e1_dp*t500*t1612 + 0.12e2_dp*t78*t168* &
    1184              :                           tnorm_drho*t91*t427 + 0.2e1_dp*t162*t163*tnorm_drhorhorho &
    1185            0 :                           - t175*t404*t553 + 0.2e1_dp*t175*t1121*t535
    1186            0 :                   t1629 = Arho*tnorm_drho1rho
    1187            0 :                   t1633 = tnorm_drho*t369
    1188            0 :                   t1646 = t361*tnorm_drho
    1189            0 :                   t1652 = t226*t369
    1190            0 :                   t1660 = t178*t1633
    1191            0 :                   t1672 = t418*tnorm_drho
    1192            0 :                   t1676 = t423*tnorm_drho
    1193              :                   t1715 = 0.2e1_dp*t999*tnorm_drho + 0.2e1_dp*t1672 + 0.2e1_dp &
    1194              :                           *t382*tnorm_drho1rho + 0.2e1_dp*t1676 + 0.2e1_dp*A*trhorho &
    1195              :                           *tnorm_drho + 0.2e1_dp*t541*tnorm_drho1rho + 0.2e1_dp*t385* &
    1196              :                           tnorm_drhorho + 0.2e1_dp*t388*tnorm_drhorho + 0.2e1_dp*t168* &
    1197              :                           tnorm_drhorhorho + 0.8e1_dp*t1187*t517 + 0.24e2_dp*t84* &
    1198              :                           t1672 + 0.8e1_dp*t417*t1629 + 0.8e1_dp*t417*t1528 + &
    1199              :                           0.24e2_dp*t84*t1676 + 0.24e2_dp*t1196*t548*t1rho + &
    1200              :                           0.12e2_dp*t426*tnorm_drho1rho*trho + 0.12e2_dp*t426* &
    1201              :                           tnorm_drho*trhorho + 0.8e1_dp*t417*t1482 + 0.12e2_dp*t426* &
    1202            0 :                           tnorm_drhorho*t1rho + 0.4e1_dp*t183*tnorm_drhorhorho
    1203              :                   t1722 = 0.2e1_dp*t217*t1629*t91 - 0.2e1_dp*t162*t397* &
    1204              :                           t1633 - t175*t397*t1565 - 0.2e1_dp*t162*t178*t226* &
    1205              :                           trhorho + 0.2e1_dp*t78*trhorho*t214 + 0.2e1_dp*t500*t1646 &
    1206              :                           - 0.2e1_dp*t217*t1521*t525 + 0.2e1_dp*t175*t1109*t1652 - &
    1207              :                           0.2e1_dp*t162*t178*tnorm_drho1rho*t186 - 0.2e1_dp*t500* &
    1208              :                           t1660 + 0.4e1_dp*t980*t1608*t400 + 0.2e1_dp*t175*t409* &
    1209              :                           t226*t432 - t175*t178*t1715 - 0.2e1_dp*t217*t218*t177* &
    1210            0 :                           t432
    1211              :                   t1758 = 0.2e1_dp*t354*t214 + 0.2e1_dp*t162*t1646 - &
    1212              :                           0.2e1_dp*t162*t1660 + 0.2e1_dp*t162*t1612 + 0.6e1_dp*t175 &
    1213              :                           *t218*t1573 + 0.2e1_dp*t217*t1511*t91 + 0.2e1_dp*t217* &
    1214              :                           t1521*t91 - 0.2e1_dp*t217*t218*t1532 - 0.2e1_dp*t162* &
    1215              :                           t178*t1515 - t175*t404*t226 + 0.2e1_dp*t175*t409*t1652 - &
    1216            0 :                           t175*t178*t1565
    1217            0 :                   t1759 = gamma_var*t1758
    1218            0 :                   t1761 = t1295*t189
    1219            0 :                   snorm_drho1rho = snorm_drhorho
    1220            0 :                   t1797 = snorm_drhorho*t103
    1221              :                   Fxnorm_drho1rho = -0.8e1_dp*t471*t566*s1rho + 0.2e1_dp* &
    1222            0 :                                     t204*s1rho*snorm_drho + 0.2e1_dp*t204*s*snorm_drho1rho
    1223              : 
    1224              :                   e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) + &
    1225              :                                         scale_ex*(ex_unif1rho*Fxnorm_drho + &
    1226              :                                                   ex_unif*Fxnorm_drho1rho + ex_unifrho*Fxnorm_drho + t487* &
    1227              :                                                   Fxnorm_drho + t208*Fxnorm_drho1rho + ex_unif*Fxnorm_drhorho + &
    1228              :                                                   t491*Fxnorm_drhorho + t108*(0.48e2_dp*t1385*snorm_drho* &
    1229              :                                                                            t1387*t476 - 0.24e2_dp*t1393*t566*t476 - 0.8e1_dp*t471* &
    1230              :                                                                            snorm_drho1rho*t103*srho - 0.8e1_dp*t471*t566*srhorho + &
    1231              :                                                                             0.2e1_dp*t204*srhorho*snorm_drho + 0.2e1_dp*t204*srho* &
    1232              :                                                                        snorm_drho1rho - 0.8e1_dp*t471*t1797*s1rho + 0.2e1_dp*t204* &
    1233              :                                                                              s1rho*snorm_drhorho + 0.2e1_dp*t204*s*(t456*t6*t458 + &
    1234              :                                                                           t196*t115*kfrho - t562*kfrhorho/0.2e1_dp + t98*t246))) + &
    1235              :                                         scale_ec*(t1759*t191 - t230*t449 + Hnorm_drhorho + my_rho*( &
    1236              :                                                   gamma_var*(t1525 + t1584 + t1628 + t1722)*t191 - t557*t449 - &
    1237            0 :                                                   t1759*t559 + 0.2e1_dp*t230*t1761*t448 - t230*t439*t435))
    1238              : 
    1239            0 :                   t1838 = t1504*t535
    1240            0 :                   t1851 = Arho*t581
    1241              :                   t1878 = 0.4e1_dp*t175*t409*t226*t553 - 0.2e1_dp*t162* &
    1242              :                           t178*t605*trho - 0.4e1_dp*t217*t218*t177*t553 + 0.8e1_dp &
    1243              :                           *t1545*t1838 + 0.2e1_dp*t78*t581*t171*t91 - 0.4e1_dp* &
    1244              :                           t162*t397*t590 - 0.10e2_dp*t175*t586*t525 - t175*t178* &
    1245              :                           (0.2e1_dp*t1851 + 0.4e1_dp*t521*tnorm_drho + 0.24e2_dp*t84* &
    1246              :                            t1851 + 0.24e2_dp*t1196*t581*trho + 0.24e2_dp*t426* &
    1247              :                            tnorm_drhorho*tnorm_drho) - 0.12e2_dp*t1492*t1493*t529 + &
    1248              :                           0.8e1_dp*t980*t1838 - 0.4e1_dp*t162*t178*tnorm_drhorho* &
    1249            0 :                           t226 + 0.20e2_dp*t162*t586*t513
    1250            0 :                   t1889 = t78*t581
    1251            0 :                   t1907 = t85*t1062
    1252              :                   t1922 = -0.4e1_dp*t500*t591 + 0.2e1_dp*t175*t409*t605* &
    1253              :                           t186 + 0.4e1_dp*t162*t409*t598*trho - 0.2e1_dp*t1889* &
    1254              :                           t187 + 0.2e1_dp*t175*t1109*t598 + 0.20e2_dp*t175*t218* &
    1255              :                           t91*tnorm_drhorho + 0.10e2_dp*t175*t1851*t91 - 0.4e1_dp* &
    1256              :                           t217*t517*t594 - t175*t397*t605 - 0.6e1_dp*t175*t1907* &
    1257              :                           t598*t186 - 0.4e1_dp*t162*t178*tnorm_drho*t553 + 0.4e1_dp &
    1258            0 :                           *t78*tnorm_drhorho*t214 - 0.4e1_dp*t217*t521*t594
    1259            0 :                   t1927 = t439*t229
    1260            0 :                   t1933 = t614*t1387
    1261            0 :                   t1937 = t614*t103
    1262              : 
    1263              :                   e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) + &
    1264              :                                           scale_ex*(ex_unif* &
    1265              :                                                     Fxnorm_drhonorm_drho + t208*Fxnorm_drhonorm_drho + t108*( &
    1266              :                                                     0.48e2_dp*t1385*t1933*srho - 0.24e2_dp*t1393*t1937*srho &
    1267              :                                                     - 0.16e2_dp*t471*t1797*snorm_drho + 0.4e1_dp*t204* &
    1268              :                                                     snorm_drhorho*snorm_drho)) + scale_ec*(Hnorm_drhonorm_drho + my_rho &
    1269              :                                                                           *(gamma_var*(t1878 + t1922)*t191 - t609*t559 - 0.2e1_dp* &
    1270            0 :                                                                                              t557*t1927 + 0.2e1_dp*t612*t1761))
    1271              : 
    1272            0 :                   t2norm_drho = tnorm_drho
    1273            0 :                   t1952 = t226*t2norm_drho
    1274            0 :                   t1964 = 0.2e1_dp*t168*t2norm_drho + 0.4e1_dp*t183*t2norm_drho
    1275            0 :                   t1965 = t178*t1964
    1276            0 :                   t1968 = t226*t1964
    1277            0 :                   t1969 = t1504*t1968
    1278            0 :                   t1972 = A*t2norm_drho
    1279              :                   t1990 = 0.2e1_dp*t1972*tnorm_drho + 0.12e2_dp*t426* &
    1280            0 :                           tnorm_drho*t2norm_drho
    1281            0 :                   t2020 = t177*t1964
    1282            0 :                   t2024 = t91*t2norm_drho
    1283            0 :                   t2028 = t78*t2norm_drho
    1284              :                   t2031 = -0.20e2_dp*t1492*t1493*t1952 + 0.4e1_dp*t162* &
    1285              :                           t409*t598*t2norm_drho - 0.2e1_dp*t1889*t1965 + 0.8e1_dp* &
    1286              :                           t1545*t1969 + 0.4e1_dp*t217*t1972*t408*t598 - 0.6e1_dp* &
    1287              :                           t175*t1907*t598*t1964 - 0.2e1_dp*t217*t1972*t177*t605 &
    1288              :                           - 0.4e1_dp*t217*t218*t177*t1990 + 0.4e1_dp*t175*t409* &
    1289              :                           t1990*t226 - 0.24e2_dp*t78*t182*t85*t177*t87*t581* &
    1290              :                           t2norm_drho + 0.2e1_dp*t175*t409*t605*t1964 + 0.8e1_dp* &
    1291              :                           t980*t1969 - 0.2e1_dp*t162*t178*t605*t2norm_drho - &
    1292              :                           0.4e1_dp*t162*t178*t1990*tnorm_drho - 0.10e2_dp*t175* &
    1293              :                           t586*t2020 + 0.24e2_dp*t162*t586*t2024 - 0.4e1_dp*t2028* &
    1294            0 :                           t591
    1295              :                   t2041 = 0.2e1_dp*t162*t163*t2norm_drho + 0.2e1_dp*t217* &
    1296            0 :                           t1972*t91 - t175*t1965
    1297            0 :                   s2norm_drho = snorm_drho
    1298              : 
    1299              :                   e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) + &
    1300              :                                             scale_ex*t108*(0.48e2_dp* &
    1301              :                                                            t1385*t1933*s2norm_drho - 0.24e2_dp*t1393*t1937* &
    1302              :                                                            s2norm_drho) + scale_ec*my_rho*(gamma_var*t2031*t191 - t609* &
    1303              :                                                                             t439*t2041 - 0.2e1_dp*gamma_var*(0.2e1_dp*t2028*t214 + &
    1304              :                                                                                    0.10e2_dp*t175*t218*t2024 - 0.2e1_dp*t162*t178* &
    1305              :                                                                            tnorm_drho*t1964 - 0.2e1_dp*t217*t218*t2020 - 0.2e1_dp* &
    1306              :                                                                             t162*t178*t1952 - 0.2e1_dp*t217*t1972*t594 + 0.2e1_dp* &
    1307              :                                                                           t175*t409*t1968 - t175*t178*t1990)*t1927 + 0.2e1_dp*t612 &
    1308            0 :                                                                                            *t1295*t2041)
    1309              :                END IF
    1310              :             END IF
    1311              :          END DO
    1312              : !$OMP END DO
    1313              :       END SELECT
    1314              : 
    1315        89699 :    END SUBROUTINE pbe_lda_calc
    1316              : 
    1317              : ! **************************************************************************************************
    1318              : !> \brief evaluates the becke 88 exchange functional for lsd
    1319              : !> \param rho_set ...
    1320              : !> \param deriv_set ...
    1321              : !> \param grad_deriv ...
    1322              : !> \param pbe_params ...
    1323              : !> \author fawzi
    1324              : ! **************************************************************************************************
    1325        96315 :    SUBROUTINE pbe_lsd_eval(rho_set, deriv_set, grad_deriv, pbe_params)
    1326              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
    1327              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
    1328              :       INTEGER, INTENT(in)                                :: grad_deriv
    1329              :       TYPE(section_vals_type), POINTER                   :: pbe_params
    1330              : 
    1331              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pbe_lsd_eval'
    1332              : 
    1333              :       INTEGER                                            :: handle, npoints, param
    1334              :       INTEGER, DIMENSION(2, 3)                           :: bo
    1335              :       REAL(kind=dp)                                      :: epsilon_rho, scale_ec, scale_ex
    1336        19263 :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, &
    1337        19263 :          e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndrb, e_ndrb_ndrb, e_ndrb_rb, e_ra, &
    1338        19263 :          e_ra_ra, e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, norm_drhob, rhoa, rhob
    1339              :       TYPE(xc_derivative_type), POINTER                  :: deriv
    1340              : 
    1341        19263 :       CALL timeset(routineN, handle)
    1342        19263 :       NULLIFY (deriv)
    1343              : 
    1344              :       CALL xc_rho_set_get(rho_set, &
    1345              :                           rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
    1346              :                           norm_drhob=norm_drhob, norm_drho=norm_drho, &
    1347              :                           rho_cutoff=epsilon_rho, &
    1348        19263 :                           local_bounds=bo)
    1349        19263 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
    1350              : 
    1351        19263 :       dummy => rhoa
    1352        19263 :       e_0 => dummy
    1353        19263 :       e_ra => dummy
    1354        19263 :       e_rb => dummy
    1355        19263 :       e_ndra_ra => dummy
    1356        19263 :       e_ndrb_rb => dummy
    1357        19263 :       e_ndr_ndr => dummy
    1358        19263 :       e_ndra_ndra => dummy
    1359        19263 :       e_ndrb_ndrb => dummy
    1360        19263 :       e_ndr => dummy
    1361        19263 :       e_ndra => dummy
    1362        19263 :       e_ndrb => dummy
    1363        19263 :       e_ra_ra => dummy
    1364        19263 :       e_ra_rb => dummy
    1365        19263 :       e_rb_rb => dummy
    1366        19263 :       e_ndr_ra => dummy
    1367        19263 :       e_ndr_rb => dummy
    1368              : 
    1369        19263 :       IF (grad_deriv >= 0) THEN
    1370              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
    1371        19263 :                                          allocate_deriv=.TRUE.)
    1372        19263 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
    1373              :       END IF
    1374        19263 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1375              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
    1376        16959 :                                          allocate_deriv=.TRUE.)
    1377        16959 :          CALL xc_derivative_get(deriv, deriv_data=e_ra)
    1378              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
    1379        16959 :                                          allocate_deriv=.TRUE.)
    1380        16959 :          CALL xc_derivative_get(deriv, deriv_data=e_rb)
    1381              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
    1382        16959 :                                          allocate_deriv=.TRUE.)
    1383        16959 :          CALL xc_derivative_get(deriv, deriv_data=e_ndr)
    1384              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
    1385        16959 :                                          allocate_deriv=.TRUE.)
    1386        16959 :          CALL xc_derivative_get(deriv, deriv_data=e_ndra)
    1387              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
    1388        16959 :                                          allocate_deriv=.TRUE.)
    1389        16959 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
    1390              :       END IF
    1391        19263 :       IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
    1392              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
    1393         1134 :                                          allocate_deriv=.TRUE.)
    1394         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ra_ra)
    1395              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
    1396         1134 :                                          allocate_deriv=.TRUE.)
    1397         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ra_rb)
    1398              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
    1399         1134 :                                          allocate_deriv=.TRUE.)
    1400         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_rb_rb)
    1401              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
    1402         1134 :                                          allocate_deriv=.TRUE.)
    1403         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ndr_ra)
    1404              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
    1405         1134 :                                          allocate_deriv=.TRUE.)
    1406         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ndr_rb)
    1407              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
    1408         1134 :                                          allocate_deriv=.TRUE.)
    1409         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ndra_ra)
    1410              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
    1411         1134 :                                          allocate_deriv=.TRUE.)
    1412         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrb_rb)
    1413              :          deriv => xc_dset_get_derivative(deriv_set, &
    1414         1134 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
    1415         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ndr_ndr)
    1416              :          deriv => xc_dset_get_derivative(deriv_set, &
    1417         1134 :                                          [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
    1418         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndra)
    1419              :          deriv => xc_dset_get_derivative(deriv_set, &
    1420         1134 :                                          [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
    1421         1134 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrb_ndrb)
    1422              :       END IF
    1423              :       IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
    1424              :          ! to do
    1425              :       END IF
    1426              : 
    1427        19263 :       CALL section_vals_val_get(pbe_params, "scale_c", r_val=scale_ec)
    1428        19263 :       CALL section_vals_val_get(pbe_params, "scale_x", r_val=scale_ex)
    1429        19263 :       CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
    1430              : 
    1431              : !$OMP PARALLEL DEFAULT(NONE),&
    1432              : !$OMP SHARED(rhoa,rhob,norm_drho,norm_drhoa,norm_drhob,e_0,e_ra,e_rb,e_ndra_ra),&
    1433              : !$OMP SHARED(e_ndrb_rb,e_ndr_ndr,e_ndra_ndra,e_ndrb_ndrb,e_ndr,e_ndra,e_ndrb),&
    1434              : !$OMP SHARED(e_ra_ra,e_ra_rb,e_rb_rb,e_ndr_ra,e_ndr_rb,grad_deriv,npoints),&
    1435        19263 : !$OMP SHARED(epsilon_rho,param,scale_ec,scale_ex)
    1436              :       CALL pbe_lsd_calc( &
    1437              :          rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa, &
    1438              :          norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb, &
    1439              :          e_ra_ndra=e_ndra_ra, &
    1440              :          e_rb_ndrb=e_ndrb_rb, e_ndr_ndr=e_ndr_ndr, &
    1441              :          e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, e_ndr=e_ndr, &
    1442              :          e_ndra=e_ndra, e_ndrb=e_ndrb, e_ra_ra=e_ra_ra, &
    1443              :          e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, e_ra_ndr=e_ndr_ra, &
    1444              :          e_rb_ndr=e_ndr_rb, &
    1445              :          grad_deriv=grad_deriv, npoints=npoints, &
    1446              :          epsilon_rho=epsilon_rho, &
    1447              :          param=param, scale_ec=scale_ec, scale_ex=scale_ex)
    1448              : !$OMP END PARALLEL
    1449              : 
    1450        19263 :       CALL timestop(handle)
    1451        19263 :    END SUBROUTINE pbe_lsd_eval
    1452              : 
    1453              : ! **************************************************************************************************
    1454              : !> \brief low level calculation of the pbe exchange-correlation functional for lsd
    1455              : !> \param rhoa ,rhob: alpha or beta spin density
    1456              : !> \param rhob ...
    1457              : !> \param norm_drho ...
    1458              : !> \param norm_drhoa ,norm_drhob,norm_drho: || grad rhoa |||,| grad rhoa ||,
    1459              : !>        || grad (rhoa+rhob) ||
    1460              : !> \param norm_drhob ...
    1461              : !> \param e_0 adds to it the local value of the functional
    1462              : !> \param e_ra derivative of the functional wrt. to the variables
    1463              : !>        named where the * is
    1464              : !> \param e_rb ...
    1465              : !> \param e_ra_ndra ...
    1466              : !> \param e_rb_ndrb ...
    1467              : !> \param e_ndr_ndr ...
    1468              : !> \param e_ndra_ndra ...
    1469              : !> \param e_ndrb_ndrb ...
    1470              : !> \param e_ndr ...
    1471              : !> \param e_ndra ...
    1472              : !> \param e_ndrb ...
    1473              : !> \param e_ra_ra ...
    1474              : !> \param e_ra_rb ...
    1475              : !> \param e_rb_rb ...
    1476              : !> \param e_ra_ndr ...
    1477              : !> \param e_rb_ndr ...
    1478              : !> \param grad_deriv ...
    1479              : !> \param npoints ...
    1480              : !> \param epsilon_rho ...
    1481              : !> \param param ...
    1482              : !> \param scale_ec ...
    1483              : !> \param scale_ex ...
    1484              : !> \author fawzi
    1485              : ! **************************************************************************************************
    1486        19263 :    SUBROUTINE pbe_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, &
    1487              :                            e_0, e_ra, e_rb, e_ra_ndra, e_rb_ndrb, e_ndr_ndr, &
    1488              :                            e_ndra_ndra, e_ndrb_ndrb, e_ndr, &
    1489              :                            e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ra_ndr, e_rb_ndr, &
    1490              :                            grad_deriv, npoints, epsilon_rho, param, scale_ec, scale_ex)
    1491              :       REAL(kind=dp), DIMENSION(*), INTENT(in)            :: rhoa, rhob, norm_drho, norm_drhoa, &
    1492              :                                                             norm_drhob
    1493              :       REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ra_ndra, e_rb_ndrb, &
    1494              :          e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb, e_ndr, e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, &
    1495              :          e_ra_ndr, e_rb_ndr
    1496              :       INTEGER, INTENT(in)                                :: grad_deriv, npoints
    1497              :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho
    1498              :       INTEGER, INTENT(in)                                :: param
    1499              :       REAL(kind=dp), INTENT(in)                          :: scale_ec, scale_ex
    1500              : 
    1501              :       INTEGER                                            :: ii
    1502              :       REAL(kind=dp) :: A, A1rhoa, A1rhob, A_1, A_2, A_3, alpha_1_1, alpha_1_2, alpha_1_3, alpha_c, &
    1503              :          alpha_c1rhoa, alpha_c1rhob, alpha_crhoa, alpha_crhob, Arhoa, Arhoarhoa, Arhoarhob, Arhob, &
    1504              :          Arhobrhob, beta, beta_1_1, beta_1_2, beta_1_3, beta_2_1, beta_2_2, beta_2_3, beta_3_1, &
    1505              :          beta_3_2, beta_3_3, beta_4_1, beta_4_2, beta_4_3, chi, chirhoa, chirhoarhoa, chirhoarhob, &
    1506              :          chirhob, chirhobrhob, e_c_u_0, e_c_u_01rhoa, e_c_u_01rhob, e_c_u_0rhoa, e_c_u_0rhoarhoa, &
    1507              :          e_c_u_0rhoarhob, e_c_u_0rhob, e_c_u_0rhobrhob, e_c_u_1rhoa, e_c_u_1rhob, epsilon_c_unif, &
    1508              :          epsilon_c_unif1rhoa, epsilon_c_unif1rhob
    1509              :       REAL(kind=dp) :: epsilon_c_unifrhoa, epsilon_c_unifrhoarhoa, epsilon_c_unifrhoarhob, &
    1510              :          epsilon_c_unifrhob, epsilon_c_unifrhobrhob, epsilon_cGGA, epsilon_cGGArhoa, &
    1511              :          epsilon_cGGArhob, ex_unif_a, ex_unif_a1rhoa, ex_unif_arhoa, ex_unif_b, ex_unif_b1rhob, &
    1512              :          ex_unif_brhob, f, f1rhoa, f1rhob, frhoa, frhoarhoa, frhoarhob, frhob, frhobrhob, Fx_a, &
    1513              :          Fx_a1rhoa, Fx_anorm_drhoa, Fx_arhoa, Fx_b, Fx_b1rhob, Fx_bnorm_drhob, Fx_brhob, &
    1514              :          gamma_var, Hnorm_drho, k_frhoa, k_frhoarhoa, k_frhoarhob, k_frhob, k_s, k_s1rhoa, &
    1515              :          k_s1rhob, k_srhoa, k_srhob, kappa, kf_a, kf_arhoa, kf_arhoarhoa, kf_b, kf_brhob, &
    1516              :          kf_brhobrhob, mu
    1517              :       REAL(kind=dp) :: my_norm_drho, my_norm_drhoa, my_norm_drhob, my_rho, my_rhoa, my_rhob, p_1, &
    1518              :          p_2, p_3, phi, phi1rhoa, phi1rhob, phirhoa, phirhoarhoa, phirhoarhob, phirhob, &
    1519              :          phirhobrhob, rs, rsrhoa, rsrhoarhoa, rsrhoarhob, rsrhob, rsrhobrhob, s_a, s_a1rhoa, &
    1520              :          s_anorm_drhoa, s_arhoa, s_b, s_b1rhob, s_bnorm_drhob, s_brhob, t, t1, t100, t1000, t1001, &
    1521              :          t101, t102, t103, t1031, t1033, t1038, t104, t105, t1050, t1069, t107, t1071, t108, t110, &
    1522              :          t1104, t1106, t111, t112, t113, t115, t116, t1164, t118, t119, t1193, t12, t122, t1228, &
    1523              :          t123, t1231, t124, t125, t126, t1269, t1283, t1285, t1286, t1288, t129
    1524              :       REAL(kind=dp) :: t1291, t1293, t130, t1304, t131, t1327, t133, t1330, t1342, t1346, t135, &
    1525              :          t137, t1377, t14, t140, t1411, t142, t143, t1440, t146, t1462, t1465, t147, t148, t1482, &
    1526              :          t1489, t1492, t15, t150, t1501, t1514, t1520, t1529, t153, t1550, t1552, t1555, t156, &
    1527              :          t1588, t1590, t1591, t1599, t1602, t1603, t162, t163, t1630, t1632, t1635, t1638, t164, &
    1528              :          t1640, t165, t167, t1674, t1677, t1680, t1698, t171, t1711, t1712, t1713, t1739, t1741, &
    1529              :          t1743, t1748, t175, t176, t1765, t1766, t1767, t177, t178, t179, t18, t1801, t1829, t183, &
    1530              :          t1830, t1831, t1865, t187, t1871, t1876, t1888, t190, t1901, t191
    1531              :       REAL(kind=dp) :: t192, t1922, t194, t1949, t198, t199, t1rhoa, t1rhob, t2, t20, t200, t201, &
    1532              :          t205, t21, t211, t212, t213, t215, t219, t22, t220, t221, t222, t226, t23, t232, t233, &
    1533              :          t234, t240, t242, t244, t245, t246, t248, t249, t250, t252, t254, t256, t257, t259, t262, &
    1534              :          t266, t268, t269, t27, t272, t273, t274, t277, t278, t28, t280, t281, t282, t284, t285, &
    1535              :          t286, t289, t293, t294, t297, t298, t299, t3, t302, t303, t305, t306, t310, t311, t312, &
    1536              :          t313, t314, t317, t318, t32, t321, t324, t325, t326, t329, t335, t336, t337, t34, t340, &
    1537              :          t341, t344, t346, t350, t368, t38, t382, t39, t396, t4, t40
    1538              :       REAL(kind=dp) :: t403, t405, t407, t409, t41, t410, t411, t413, t415, t417, t427, t429, &
    1539              :          t432, t436, t439, t442, t444, t445, t45, t453, t456, t457, t46, t460, t466, t467, t468, &
    1540              :          t471, t472, t475, t477, t481, t488, t493, t494, t5, t50, t502, t505, t506, t518, t519, &
    1541              :          t52, t523, t525, t536, t538, t543, t544, t548, t549, t550, t556, t56, t561, t563, t564, &
    1542              :          t57, t576, t578, t579, t58, t580, t588, t59, t590, t595, t596, t6, t600, t606, t611, &
    1543              :          t624, t626, t627, t628, t63, t636, t638, t64, t643, t644, t648, t654, t659, t66, t672, &
    1544              :          t674, t675, t676, t681, t682, t687, t69, t7, t70, t705, t708, t71, t711
    1545              :       REAL(kind=dp) :: t72, t726, t73, t733, t736, t74, t745, t75, t750, t755, t763, t767, t768, &
    1546              :          t77, t775, t776, t779, t78, t785, t789, t79, t795, t798, t8, t80, t801, t812, t82, t820, &
    1547              :          t821, t822, t823, t825, t828, t839, t84, t840, t85, t851, t858, t865, t867, t868, t87, &
    1548              :          t876, t879, t88, t880, t9, t90, t904, t908, t909, t91, t911, t914, t917, t919, t92, t924, &
    1549              :          t93, t936, t94, t944, t95, t953, t959, t96, t962, t965, t966, t967, t97, t98, t985, t998, &
    1550              :          t999, tnorm_drho, trhoa, trhoanorm_drho, trhoarhoa, trhoarhob, trhob, trhobnorm_drho, &
    1551              :          trhobrhob
    1552              : 
    1553              : ! Parametrization of PBE
    1554              : 
    1555        19263 :       SELECT CASE (param)
    1556              :          ! Original PBE
    1557              :       CASE (xc_pbe_orig)
    1558              :          kappa = 0.804e0_dp
    1559              :          beta = 0.66725e-1_dp
    1560           20 :          mu = beta*pi**2/0.3e1_dp
    1561              :          ! Revised PBE (revPBE)
    1562              :       CASE (xc_pbe_rev)
    1563           20 :          kappa = 0.1245e1_dp
    1564           20 :          beta = 0.66725e-1_dp
    1565           20 :          mu = beta*pi**2/0.3e1_dp
    1566              :          ! PBE for solids (PBEsol)
    1567              :       CASE (xc_pbe_sol)
    1568          142 :          kappa = 0.804e0_dp
    1569          142 :          beta = 0.46e-1_dp
    1570          142 :          mu = 0.1e2_dp/0.81e2_dp
    1571              :       CASE default
    1572        19263 :          CPABORT("Unsupported parametrization")
    1573              :       END SELECT
    1574              : 
    1575        19263 :       gamma_var = (0.1e1_dp - LOG(0.2e1_dp))/pi**2
    1576        19263 :       p_1 = 0.10e1_dp
    1577        19263 :       A_1 = 0.31091e-1_dp
    1578        19263 :       alpha_1_1 = 0.21370e0_dp
    1579        19263 :       beta_1_1 = 0.75957e1_dp
    1580        19263 :       beta_2_1 = 0.35876e1_dp
    1581        19263 :       beta_3_1 = 0.16382e1_dp
    1582        19263 :       beta_4_1 = 0.49294e0_dp
    1583        19263 :       p_2 = 0.10e1_dp
    1584        19263 :       A_2 = 0.15545e-1_dp
    1585        19263 :       alpha_1_2 = 0.20548e0_dp
    1586        19263 :       beta_1_2 = 0.141189e2_dp
    1587        19263 :       beta_2_2 = 0.61977e1_dp
    1588        19263 :       beta_3_2 = 0.33662e1_dp
    1589        19263 :       beta_4_2 = 0.62517e0_dp
    1590        19263 :       p_3 = 0.10e1_dp
    1591        19263 :       A_3 = 0.16887e-1_dp
    1592        19263 :       alpha_1_3 = 0.11125e0_dp
    1593        19263 :       beta_1_3 = 0.10357e2_dp
    1594        19263 :       beta_2_3 = 0.36231e1_dp
    1595        19263 :       beta_3_3 = 0.88026e0_dp
    1596        19263 :       beta_4_3 = 0.49671e0_dp
    1597        19263 :       t3 = 3**(0.1e1_dp/0.3e1_dp)
    1598        19263 :       t4 = 4**(0.1e1_dp/0.3e1_dp)
    1599        19263 :       t5 = t4**2
    1600        19263 :       t6 = t3*t5
    1601        19263 :       t7 = 0.1e1_dp/pi
    1602        19263 :       t90 = pi**2
    1603              : 
    1604              :       SELECT CASE (grad_deriv)
    1605              :       CASE default
    1606        19263 : !$OMP DO
    1607              :          DO ii = 1, npoints
    1608    624328540 :             my_rhoa = MAX(rhoa(ii), 0.0_dp)
    1609    624328540 :             my_rhob = MAX(rhob(ii), 0.0_dp)
    1610    624328540 :             my_rho = my_rhoa + my_rhob
    1611    624328540 :             IF (my_rho > epsilon_rho) THEN
    1612    580030897 :                my_rhoa = MAX(EPSILON(0.0_dp)*1.e4_dp, my_rhoa)
    1613    580030897 :                my_rhob = MAX(EPSILON(0.0_dp)*1.e4_dp, my_rhob)
    1614    580030897 :                my_rho = my_rhoa + my_rhob
    1615    580030897 :                my_norm_drho = norm_drho(ii)
    1616    580030897 :                my_norm_drhoa = norm_drhoa(ii)
    1617    580030897 :                my_norm_drhob = norm_drhob(ii)
    1618              : 
    1619    580030897 :                t1 = my_rhoa - my_rhob
    1620    580030897 :                t2 = 0.1e1_dp/my_rho
    1621    580030897 :                chi = t1*t2
    1622    580030897 :                t8 = t7*t2
    1623    580030897 :                t9 = t8**(0.1e1_dp/0.3e1_dp)
    1624    580030897 :                rs = t6*t9/0.4e1_dp
    1625    580030897 :                t12 = 0.1e1_dp + alpha_1_1*rs
    1626    580030897 :                t14 = 0.1e1_dp/A_1
    1627    580030897 :                t15 = SQRT(rs)
    1628    580030897 :                t18 = t15*rs
    1629    580030897 :                t20 = p_1 + 0.1e1_dp
    1630    580030897 :                t21 = rs**t20
    1631    580030897 :                t22 = beta_4_1*t21
    1632    580030897 :                t23 = beta_1_1*t15 + beta_2_1*rs + beta_3_1*t18 + t22
    1633    580030897 :                t27 = 0.1e1_dp + t14/t23/0.2e1_dp
    1634    580030897 :                t28 = LOG(t27)
    1635    580030897 :                e_c_u_0 = -0.2e1_dp*A_1*t12*t28
    1636    580030897 :                t32 = 0.1e1_dp + alpha_1_2*rs
    1637    580030897 :                t34 = 0.1e1_dp/A_2
    1638    580030897 :                t38 = p_2 + 0.1e1_dp
    1639    580030897 :                t39 = rs**t38
    1640    580030897 :                t40 = beta_4_2*t39
    1641    580030897 :                t41 = beta_1_2*t15 + beta_2_2*rs + beta_3_2*t18 + t40
    1642    580030897 :                t45 = 0.1e1_dp + t34/t41/0.2e1_dp
    1643    580030897 :                t46 = LOG(t45)
    1644    580030897 :                t50 = 0.1e1_dp + alpha_1_3*rs
    1645    580030897 :                t52 = 0.1e1_dp/A_3
    1646    580030897 :                t56 = p_3 + 0.1e1_dp
    1647    580030897 :                t57 = rs**t56
    1648    580030897 :                t58 = beta_4_3*t57
    1649    580030897 :                t59 = beta_1_3*t15 + beta_2_3*rs + beta_3_3*t18 + t58
    1650    580030897 :                t63 = 0.1e1_dp + t52/t59/0.2e1_dp
    1651    580030897 :                t64 = LOG(t63)
    1652    580030897 :                alpha_c = 0.2e1_dp*A_3*t50*t64
    1653    580030897 :                t66 = 2**(0.1e1_dp/0.3e1_dp)
    1654    580030897 :                t69 = 1/(2*t66 - 2)
    1655    580030897 :                t70 = 0.1e1_dp + chi
    1656    580030897 :                t71 = t70**(0.1e1_dp/0.3e1_dp)
    1657    580030897 :                t72 = t71*t70
    1658    580030897 :                t73 = 0.1e1_dp - chi
    1659    580030897 :                t74 = t73**(0.1e1_dp/0.3e1_dp)
    1660    580030897 :                t75 = t74*t73
    1661    580030897 :                f = (t72 + t75 - 0.2e1_dp)*t69
    1662    580030897 :                t77 = alpha_c*f
    1663    580030897 :                t78 = 0.9e1_dp/0.8e1_dp/t69
    1664    580030897 :                t79 = chi**2
    1665    580030897 :                t80 = t79**2
    1666    580030897 :                t82 = t78*(0.1e1_dp - t80)
    1667    580030897 :                t84 = -0.2e1_dp*A_2*t32*t46 - e_c_u_0
    1668    580030897 :                t85 = t84*f
    1669    580030897 :                epsilon_c_unif = e_c_u_0 + t77*t82 + t85*t80
    1670    580030897 :                t87 = t71**2
    1671    580030897 :                t88 = t74**2
    1672    580030897 :                phi = t87/0.2e1_dp + t88/0.2e1_dp
    1673    580030897 :                t91 = t90*my_rho
    1674    580030897 :                t92 = t91**(0.1e1_dp/0.3e1_dp)
    1675    580030897 :                t93 = t3*t92*t7
    1676    580030897 :                t94 = SQRT(t93)
    1677    580030897 :                k_s = 0.2e1_dp*t94
    1678    580030897 :                t95 = 0.1e1_dp/phi
    1679    580030897 :                t96 = my_norm_drho*t95
    1680    580030897 :                t97 = 0.1e1_dp/k_s
    1681    580030897 :                t98 = t97*t2
    1682    580030897 :                t = t96*t98/0.2e1_dp
    1683    580030897 :                t100 = 0.1e1_dp/gamma_var
    1684    580030897 :                t101 = beta*t100
    1685    580030897 :                t102 = epsilon_c_unif*t100
    1686    580030897 :                t103 = phi**2
    1687    580030897 :                t104 = t103*phi
    1688    580030897 :                t105 = 0.1e1_dp/t104
    1689    580030897 :                t107 = EXP(-t102*t105)
    1690    580030897 :                t108 = t107 - 0.1e1_dp
    1691    580030897 :                A = t101/t108
    1692    580030897 :                t110 = gamma_var*t104
    1693    580030897 :                t111 = t**2
    1694    580030897 :                t112 = A*t111
    1695    580030897 :                t113 = 0.1e1_dp + t112
    1696    580030897 :                t115 = A**2
    1697    580030897 :                t116 = t111**2
    1698    580030897 :                t118 = 0.1e1_dp + t112 + t115*t116
    1699    580030897 :                t119 = 0.1e1_dp/t118
    1700    580030897 :                t122 = 0.1e1_dp + t101*t111*t113*t119
    1701    580030897 :                t123 = LOG(t122)
    1702    580030897 :                epsilon_cGGA = epsilon_c_unif + t110*t123
    1703    580030897 :                t124 = t3*t66
    1704    580030897 :                t125 = t90*my_rhoa
    1705    580030897 :                t126 = t125**(0.1e1_dp/0.3e1_dp)
    1706    580030897 :                kf_a = t124*t126
    1707    580030897 :                ex_unif_a = -0.3e1_dp/0.4e1_dp*t7*kf_a
    1708    580030897 :                t129 = 0.1e1_dp/kf_a
    1709    580030897 :                t130 = my_norm_drhoa*t129
    1710    580030897 :                t131 = 0.1e1_dp/my_rhoa
    1711    580030897 :                s_a = t130*t131/0.2e1_dp
    1712    580030897 :                t133 = s_a**2
    1713    580030897 :                t135 = 0.1e1_dp/kappa
    1714    580030897 :                t137 = 0.1e1_dp + mu*t133*t135
    1715    580030897 :                Fx_a = 0.1e1_dp + kappa - kappa/t137
    1716    580030897 :                t140 = my_rhoa*ex_unif_a
    1717    580030897 :                t142 = t90*my_rhob
    1718    580030897 :                t143 = t142**(0.1e1_dp/0.3e1_dp)
    1719    580030897 :                kf_b = t124*t143
    1720    580030897 :                ex_unif_b = -0.3e1_dp/0.4e1_dp*t7*kf_b
    1721    580030897 :                t146 = 0.1e1_dp/kf_b
    1722    580030897 :                t147 = my_norm_drhob*t146
    1723    580030897 :                t148 = 0.1e1_dp/my_rhob
    1724    580030897 :                s_b = t147*t148/0.2e1_dp
    1725    580030897 :                t150 = s_b**2
    1726    580030897 :                t153 = 0.1e1_dp + mu*t150*t135
    1727    580030897 :                Fx_b = 0.1e1_dp + kappa - kappa/t153
    1728    580030897 :                t156 = my_rhob*ex_unif_b
    1729              : 
    1730    580030897 :                IF (grad_deriv >= 0) THEN
    1731              :                   e_0(ii) = e_0(ii) + &
    1732              :                             scale_ex*(0.2e1_dp*t140*Fx_a + 0.2e1_dp*t156*Fx_b) &
    1733    580030897 :                             /0.2e1_dp + scale_ec*my_rho*epsilon_cGGA
    1734              :                END IF
    1735              : 
    1736    580030897 :                t162 = my_rho**2
    1737    580030897 :                t163 = 0.1e1_dp/t162
    1738    580030897 :                t164 = t1*t163
    1739    580030897 :                chirhoa = t2 - t164
    1740    580030897 :                t165 = t9**2
    1741    580030897 :                t167 = 0.1e1_dp/t165*t7
    1742    580030897 :                rsrhoa = -t6*t167*t163/0.12e2_dp
    1743    580030897 :                t171 = A_1*alpha_1_1
    1744    580030897 :                t175 = t23**2
    1745    580030897 :                t176 = 0.1e1_dp/t175
    1746    580030897 :                t177 = t12*t176
    1747    580030897 :                t178 = 0.1e1_dp/t15
    1748    580030897 :                t179 = beta_1_1*t178
    1749    580030897 :                t183 = beta_3_1*t15
    1750    580030897 :                t187 = 0.1e1_dp/rs
    1751              :                t190 = t179*rsrhoa/0.2e1_dp + beta_2_1*rsrhoa + 0.3e1_dp/ &
    1752    580030897 :                       0.2e1_dp*t183*rsrhoa + t22*t20*rsrhoa*t187
    1753    580030897 :                t191 = 0.1e1_dp/t27
    1754    580030897 :                t192 = t190*t191
    1755    580030897 :                e_c_u_0rhoa = -0.2e1_dp*t171*rsrhoa*t28 + t177*t192
    1756    580030897 :                t194 = A_2*alpha_1_2
    1757    580030897 :                t198 = t41**2
    1758    580030897 :                t199 = 0.1e1_dp/t198
    1759    580030897 :                t200 = t32*t199
    1760    580030897 :                t201 = beta_1_2*t178
    1761    580030897 :                t205 = beta_3_2*t15
    1762              :                t211 = t201*rsrhoa/0.2e1_dp + beta_2_2*rsrhoa + 0.3e1_dp/ &
    1763    580030897 :                       0.2e1_dp*t205*rsrhoa + t40*t38*rsrhoa*t187
    1764    580030897 :                t212 = 0.1e1_dp/t45
    1765    580030897 :                t213 = t211*t212
    1766    580030897 :                e_c_u_1rhoa = -0.2e1_dp*t194*rsrhoa*t46 + t200*t213
    1767    580030897 :                t215 = A_3*alpha_1_3
    1768    580030897 :                t219 = t59**2
    1769    580030897 :                t220 = 0.1e1_dp/t219
    1770    580030897 :                t221 = t50*t220
    1771    580030897 :                t222 = beta_1_3*t178
    1772    580030897 :                t226 = beta_3_3*t15
    1773              :                t232 = t222*rsrhoa/0.2e1_dp + beta_2_3*rsrhoa + 0.3e1_dp/ &
    1774    580030897 :                       0.2e1_dp*t226*rsrhoa + t58*t56*rsrhoa*t187
    1775    580030897 :                t233 = 0.1e1_dp/t63
    1776    580030897 :                t234 = t232*t233
    1777    580030897 :                alpha_crhoa = 0.2e1_dp*t215*rsrhoa*t64 - t221*t234
    1778              :                frhoa = (0.4e1_dp/0.3e1_dp*t71*chirhoa - 0.4e1_dp/0.3e1_dp &
    1779    580030897 :                         *t74*chirhoa)*t69
    1780    580030897 :                t240 = alpha_crhoa*f
    1781    580030897 :                t242 = alpha_c*frhoa
    1782    580030897 :                t244 = t79*chi
    1783    580030897 :                t245 = t78*t244
    1784    580030897 :                t246 = t245*chirhoa
    1785    580030897 :                t248 = 0.4e1_dp*t77*t246
    1786    580030897 :                t249 = e_c_u_1rhoa - e_c_u_0rhoa
    1787    580030897 :                t250 = t249*f
    1788    580030897 :                t252 = t84*frhoa
    1789    580030897 :                t254 = t244*chirhoa
    1790    580030897 :                t256 = 0.4e1_dp*t85*t254
    1791              :                epsilon_c_unifrhoa = e_c_u_0rhoa + t240*t82 + t242*t82 - t248 &
    1792    580030897 :                                     + t250*t80 + t252*t80 + t256
    1793    580030897 :                t257 = 0.1e1_dp/t71
    1794    580030897 :                t259 = 0.1e1_dp/t74
    1795    580030897 :                phirhoa = t257*chirhoa/0.3e1_dp - t259*chirhoa/0.3e1_dp
    1796    580030897 :                t262 = t92**2
    1797    580030897 :                k_frhoa = t3/t262*t90/0.3e1_dp
    1798    580030897 :                t266 = 0.1e1_dp/t94
    1799    580030897 :                k_srhoa = t266*k_frhoa*t7
    1800    580030897 :                t268 = 0.1e1_dp/t103
    1801    580030897 :                t269 = my_norm_drho*t268
    1802    580030897 :                t272 = k_s**2
    1803    580030897 :                t273 = 0.1e1_dp/t272
    1804    580030897 :                t274 = t273*t2
    1805    580030897 :                t277 = t97*t163
    1806    580030897 :                t278 = t96*t277
    1807              :                trhoa = -t269*t98*phirhoa/0.2e1_dp - t96*t274*k_srhoa/ &
    1808    580030897 :                        0.2e1_dp - t278/0.2e1_dp
    1809    580030897 :                t280 = t108**2
    1810    580030897 :                t281 = 0.1e1_dp/t280
    1811    580030897 :                t282 = epsilon_c_unifrhoa*t100
    1812    580030897 :                t284 = t103**2
    1813    580030897 :                t285 = 0.1e1_dp/t284
    1814    580030897 :                t286 = t285*phirhoa
    1815    580030897 :                t289 = -t282*t105 + 0.3e1_dp*t102*t286
    1816    580030897 :                Arhoa = -t101*t281*t289*t107
    1817    580030897 :                t293 = gamma_var*t103
    1818    580030897 :                t294 = t123*phirhoa
    1819    580030897 :                t297 = t101*t
    1820    580030897 :                t298 = t113*t119
    1821    580030897 :                t299 = t298*trhoa
    1822    580030897 :                t302 = Arhoa*t111
    1823    580030897 :                t303 = A*t
    1824    580030897 :                t305 = 0.2e1_dp*t303*trhoa
    1825    580030897 :                t306 = t302 + t305
    1826    580030897 :                t310 = t101*t111
    1827    580030897 :                t311 = t118**2
    1828    580030897 :                t312 = 0.1e1_dp/t311
    1829    580030897 :                t313 = t113*t312
    1830    580030897 :                t314 = A*t116
    1831    580030897 :                t317 = t111*t
    1832    580030897 :                t318 = t115*t317
    1833    580030897 :                t321 = t302 + t305 + 0.2e1_dp*t314*Arhoa + 0.4e1_dp*t318*trhoa
    1834              :                t324 = 0.2e1_dp*t297*t299 + t101*t111*t306*t119 - t310* &
    1835    580030897 :                       t313*t321
    1836    580030897 :                t325 = 0.1e1_dp/t122
    1837    580030897 :                t326 = t324*t325
    1838              :                epsilon_cGGArhoa = epsilon_c_unifrhoa + 0.3e1_dp*t293*t294 + &
    1839    580030897 :                                   t110*t326
    1840    580030897 :                t329 = t126**2
    1841    580030897 :                kf_arhoa = t124/t329*t90/0.3e1_dp
    1842    580030897 :                ex_unif_arhoa = -0.3e1_dp/0.4e1_dp*t7*kf_arhoa
    1843    580030897 :                t335 = kf_a**2
    1844    580030897 :                t336 = 0.1e1_dp/t335
    1845    580030897 :                t337 = my_norm_drhoa*t336
    1846    580030897 :                t340 = my_rhoa**2
    1847    580030897 :                t341 = 0.1e1_dp/t340
    1848    580030897 :                s_arhoa = -t337*t131*kf_arhoa/0.2e1_dp - t130*t341/0.2e1_dp
    1849    580030897 :                t344 = t137**2
    1850    580030897 :                t346 = 0.1e1_dp/t344*mu
    1851    580030897 :                Fx_arhoa = 0.2e1_dp*t346*s_a*s_arhoa
    1852    580030897 :                t350 = my_rhoa*ex_unif_arhoa
    1853              : 
    1854    580030897 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1855              :                   e_ra(ii) = e_ra(ii) + &
    1856              :                              scale_ex*(0.2e1_dp*ex_unif_a*Fx_a + 0.2e1_dp* &
    1857              :                                        t350*Fx_a + 0.2e1_dp*t140*Fx_arhoa)/0.2e1_dp + scale_ec*( &
    1858    476782003 :                              epsilon_cGGA + my_rho*epsilon_cGGArhoa)
    1859              :                END IF
    1860              : 
    1861    580030897 :                chirhob = -t2 - t164
    1862    580030897 :                rsrhob = rsrhoa
    1863              :                t368 = t179*rsrhob/0.2e1_dp + beta_2_1*rsrhob + 0.3e1_dp/ &
    1864    580030897 :                       0.2e1_dp*t183*rsrhob + t22*t20*rsrhob*t187
    1865    580030897 :                e_c_u_0rhob = -0.2e1_dp*t171*rsrhob*t28 + t177*t368*t191
    1866              :                t382 = t201*rsrhob/0.2e1_dp + beta_2_2*rsrhob + 0.3e1_dp/ &
    1867    580030897 :                       0.2e1_dp*t205*rsrhob + t40*t38*rsrhob*t187
    1868    580030897 :                e_c_u_1rhob = -0.2e1_dp*t194*rsrhob*t46 + t200*t382*t212
    1869              :                t396 = t222*rsrhob/0.2e1_dp + beta_2_3*rsrhob + 0.3e1_dp/ &
    1870    580030897 :                       0.2e1_dp*t226*rsrhob + t58*t56*rsrhob*t187
    1871    580030897 :                alpha_crhob = 0.2e1_dp*t215*rsrhob*t64 - t221*t396*t233
    1872              :                frhob = (0.4e1_dp/0.3e1_dp*t71*chirhob - 0.4e1_dp/0.3e1_dp &
    1873    580030897 :                         *t74*chirhob)*t69
    1874    580030897 :                t403 = alpha_crhob*f
    1875    580030897 :                t405 = alpha_c*frhob
    1876    580030897 :                t407 = t245*chirhob
    1877    580030897 :                t409 = 0.4e1_dp*t77*t407
    1878    580030897 :                t410 = e_c_u_1rhob - e_c_u_0rhob
    1879    580030897 :                t411 = t410*f
    1880    580030897 :                t413 = t84*frhob
    1881    580030897 :                t415 = t244*chirhob
    1882    580030897 :                t417 = 0.4e1_dp*t85*t415
    1883              :                epsilon_c_unifrhob = e_c_u_0rhob + t403*t82 + t405*t82 - t409 &
    1884    580030897 :                                     + t411*t80 + t413*t80 + t417
    1885    580030897 :                phirhob = t257*chirhob/0.3e1_dp - t259*chirhob/0.3e1_dp
    1886    580030897 :                k_frhob = k_frhoa
    1887    580030897 :                k_srhob = t266*k_frhob*t7
    1888              :                trhob = -t269*t98*phirhob/0.2e1_dp - t96*t274*k_srhob/ &
    1889    580030897 :                        0.2e1_dp - t278/0.2e1_dp
    1890    580030897 :                t427 = epsilon_c_unifrhob*t100
    1891    580030897 :                t429 = t285*phirhob
    1892    580030897 :                t432 = -t427*t105 + 0.3e1_dp*t102*t429
    1893    580030897 :                Arhob = -t101*t281*t432*t107
    1894    580030897 :                t436 = t123*phirhob
    1895    580030897 :                t439 = t298*trhob
    1896    580030897 :                t442 = Arhob*t111
    1897    580030897 :                t444 = 0.2e1_dp*t303*trhob
    1898    580030897 :                t445 = t442 + t444
    1899    580030897 :                t453 = t442 + t444 + 0.2e1_dp*t314*Arhob + 0.4e1_dp*t318*trhob
    1900              :                t456 = 0.2e1_dp*t297*t439 + t101*t111*t445*t119 - t310* &
    1901    580030897 :                       t313*t453
    1902    580030897 :                t457 = t456*t325
    1903              :                epsilon_cGGArhob = epsilon_c_unifrhob + 0.3e1_dp*t293*t436 + &
    1904    580030897 :                                   t110*t457
    1905    580030897 :                t460 = t143**2
    1906    580030897 :                kf_brhob = t124/t460*t90/0.3e1_dp
    1907    580030897 :                ex_unif_brhob = -0.3e1_dp/0.4e1_dp*t7*kf_brhob
    1908    580030897 :                t466 = kf_b**2
    1909    580030897 :                t467 = 0.1e1_dp/t466
    1910    580030897 :                t468 = my_norm_drhob*t467
    1911    580030897 :                t471 = my_rhob**2
    1912    580030897 :                t472 = 0.1e1_dp/t471
    1913    580030897 :                s_brhob = -t468*t148*kf_brhob/0.2e1_dp - t147*t472/0.2e1_dp
    1914    580030897 :                t475 = t153**2
    1915    580030897 :                t477 = 0.1e1_dp/t475*mu
    1916    580030897 :                Fx_brhob = 0.2e1_dp*t477*s_b*s_brhob
    1917    580030897 :                t481 = my_rhob*ex_unif_brhob
    1918              : 
    1919    580030897 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1920              :                   e_rb(ii) = e_rb(ii) + &
    1921              :                              scale_ex*(0.2e1_dp*ex_unif_b*Fx_b + 0.2e1_dp* &
    1922              :                                        t481*Fx_b + 0.2e1_dp*t156*Fx_brhob)/0.2e1_dp + scale_ec*( &
    1923    476782003 :                              epsilon_cGGA + my_rho*epsilon_cGGArhob)
    1924              :                END IF
    1925              : 
    1926    580030897 :                t488 = t95*t97
    1927    580030897 :                tnorm_drho = t488*t2/0.2e1_dp
    1928    580030897 :                t493 = t101*t317
    1929    580030897 :                t494 = A*tnorm_drho
    1930    580030897 :                t502 = 0.2e1_dp*t303*tnorm_drho + 0.4e1_dp*t318*tnorm_drho
    1931              :                t505 = 0.2e1_dp*t297*t298*tnorm_drho + 0.2e1_dp*t493* &
    1932    580030897 :                       t494*t119 - t310*t313*t502
    1933    580030897 :                t506 = t505*t325
    1934    580030897 :                Hnorm_drho = t110*t506
    1935              : 
    1936    580030897 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1937              :                   e_ndr(ii) = e_ndr(ii) + &
    1938    476782003 :                               scale_ec*my_rho*Hnorm_drho
    1939              :                END IF
    1940              : 
    1941    580030897 :                s_anorm_drhoa = t129*t131/0.2e1_dp
    1942    580030897 :                Fx_anorm_drhoa = 0.2e1_dp*t346*s_a*s_anorm_drhoa
    1943              : 
    1944    580030897 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1945              :                   e_ndra(ii) = e_ndra(ii) + &
    1946    476782003 :                                scale_ex*t140*Fx_anorm_drhoa
    1947              :                END IF
    1948              : 
    1949    580030897 :                s_bnorm_drhob = t146*t148/0.2e1_dp
    1950    580030897 :                Fx_bnorm_drhob = 0.2e1_dp*t477*s_b*s_bnorm_drhob
    1951              : 
    1952    580030897 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1953              :                   e_ndrb(ii) = e_ndrb(ii) + &
    1954    476782003 :                                scale_ex*t156*Fx_bnorm_drhob
    1955              :                END IF
    1956              : 
    1957    580030897 :                IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
    1958     18727793 :                   t518 = 0.1e1_dp/t162/my_rho
    1959     18727793 :                   t519 = t1*t518
    1960     18727793 :                   chirhoarhoa = -0.2e1_dp*t163 + 0.2e1_dp*t519
    1961     18727793 :                   t523 = 0.1e1_dp/t90
    1962     18727793 :                   t525 = t162**2
    1963              :                   rsrhoarhoa = -t6/t165/t8*t523/t525/0.18e2_dp + &
    1964     18727793 :                                t6*t167*t518/0.6e1_dp
    1965     18727793 :                   t536 = alpha_1_1*rsrhoa
    1966     18727793 :                   t538 = t176*t190*t191
    1967     18727793 :                   t543 = t12/t175/t23
    1968     18727793 :                   t544 = t190**2
    1969     18727793 :                   t548 = 0.1e1_dp/t18
    1970     18727793 :                   t549 = beta_1_1*t548
    1971     18727793 :                   t550 = rsrhoa**2
    1972     18727793 :                   t556 = beta_3_1*t178
    1973     18727793 :                   t561 = t20**2
    1974     18727793 :                   t563 = rs**2
    1975     18727793 :                   t564 = 0.1e1_dp/t563
    1976     18727793 :                   t576 = t175**2
    1977     18727793 :                   t578 = t12/t576
    1978     18727793 :                   t579 = t27**2
    1979     18727793 :                   t580 = 0.1e1_dp/t579
    1980              :                   e_c_u_0rhoarhoa = -0.2e1_dp*t171*rsrhoarhoa*t28 + 0.2e1_dp* &
    1981              :                                     t536*t538 - 0.2e1_dp*t543*t544*t191 + t177*(-t549*t550 &
    1982              :                                                                       /0.4e1_dp + t179*rsrhoarhoa/0.2e1_dp + beta_2_1*rsrhoarhoa + &
    1983              :                                                                              0.3e1_dp/0.4e1_dp*t556*t550 + 0.3e1_dp/0.2e1_dp*t183* &
    1984              :                                                                              rsrhoarhoa + t22*t561*t550*t564 + t22*t20*rsrhoarhoa* &
    1985              :                                                                               t187 - t22*t20*t550*t564)*t191 + t578*t544*t580*t14/ &
    1986     18727793 :                                     0.2e1_dp
    1987     18727793 :                   e_c_u_01rhoa = e_c_u_0rhoa
    1988     18727793 :                   t588 = alpha_1_2*rsrhoa
    1989     18727793 :                   t590 = t199*t211*t212
    1990     18727793 :                   t595 = t32/t198/t41
    1991     18727793 :                   t596 = t211**2
    1992     18727793 :                   t600 = beta_1_2*t548
    1993     18727793 :                   t606 = beta_3_2*t178
    1994     18727793 :                   t611 = t38**2
    1995     18727793 :                   t624 = t198**2
    1996     18727793 :                   t626 = t32/t624
    1997     18727793 :                   t627 = t45**2
    1998     18727793 :                   t628 = 0.1e1_dp/t627
    1999     18727793 :                   t636 = alpha_1_3*rsrhoa
    2000     18727793 :                   t638 = t220*t232*t233
    2001     18727793 :                   t643 = t50/t219/t59
    2002     18727793 :                   t644 = t232**2
    2003     18727793 :                   t648 = beta_1_3*t548
    2004     18727793 :                   t654 = beta_3_3*t178
    2005     18727793 :                   t659 = t56**2
    2006     18727793 :                   t672 = t219**2
    2007     18727793 :                   t674 = t50/t672
    2008     18727793 :                   t675 = t63**2
    2009     18727793 :                   t676 = 0.1e1_dp/t675
    2010     18727793 :                   alpha_c1rhoa = alpha_crhoa
    2011     18727793 :                   t681 = 0.1e1_dp/t87
    2012     18727793 :                   t682 = chirhoa**2
    2013     18727793 :                   t687 = 0.1e1_dp/t88
    2014              :                   frhoarhoa = (0.4e1_dp/0.9e1_dp*t681*t682 + 0.4e1_dp/ &
    2015              :                                0.3e1_dp*t71*chirhoarhoa + 0.4e1_dp/0.9e1_dp*t687*t682 - &
    2016     18727793 :                                0.4e1_dp/0.3e1_dp*t74*chirhoarhoa)*t69
    2017     18727793 :                   f1rhoa = frhoa
    2018     18727793 :                   t705 = alpha_c1rhoa*f
    2019     18727793 :                   t708 = alpha_c*f1rhoa
    2020     18727793 :                   t711 = t78*t79
    2021     18727793 :                   t726 = e_c_u_1rhoa - e_c_u_01rhoa
    2022     18727793 :                   t733 = t726*f
    2023     18727793 :                   t736 = t84*f1rhoa
    2024              :                   t745 = -0.4e1_dp*t77*t245*chirhoarhoa + (-0.2e1_dp*t194* &
    2025              :                                                            rsrhoarhoa*t46 + 0.2e1_dp*t588*t590 - 0.2e1_dp*t595*t596* &
    2026              :                                                            t212 + t200*(-t600*t550/0.4e1_dp + t201*rsrhoarhoa/ &
    2027              :                                                                       0.2e1_dp + beta_2_2*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t606*t550 &
    2028              :                                                                         + 0.3e1_dp/0.2e1_dp*t205*rsrhoarhoa + t40*t611*t550* &
    2029              :                                                                         t564 + t40*t38*rsrhoarhoa*t187 - t40*t38*t550*t564)* &
    2030              :                                                            t212 + t626*t596*t628*t34/0.2e1_dp - e_c_u_0rhoarhoa)*f* &
    2031              :                          t80 + t249*f1rhoa*t80 + 0.4e1_dp*t250*t254 + t726*frhoa* &
    2032              :                          t80 + t84*frhoarhoa*t80 + 0.4e1_dp*t252*t254 + 0.4e1_dp* &
    2033              :                          t733*t254 + 0.4e1_dp*t736*t254 + 0.12e2_dp*t85*t79*t682 &
    2034     18727793 :                          + 0.4e1_dp*t85*t244*chirhoarhoa
    2035              :                   epsilon_c_unifrhoarhoa = e_c_u_0rhoarhoa + (0.2e1_dp*t215* &
    2036              :                                                               rsrhoarhoa*t64 - 0.2e1_dp*t636*t638 + 0.2e1_dp*t643*t644* &
    2037              :                                                               t233 - t221*(-t648*t550/0.4e1_dp + t222*rsrhoarhoa/ &
    2038              :                                                                       0.2e1_dp + beta_2_3*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t654*t550 &
    2039              :                                                                            + 0.3e1_dp/0.2e1_dp*t226*rsrhoarhoa + t58*t659*t550* &
    2040              :                                                                            t564 + t58*t56*rsrhoarhoa*t187 - t58*t56*t550*t564)* &
    2041              :                                                               t233 - t674*t644*t676*t52/0.2e1_dp)*f*t82 + alpha_crhoa &
    2042              :                                            *f1rhoa*t82 - 0.4e1_dp*t240*t246 + alpha_c1rhoa*frhoa*t82 &
    2043              :                                            + alpha_c*frhoarhoa*t82 - 0.4e1_dp*t242*t246 - 0.4e1_dp* &
    2044              :                                            t705*t246 - 0.4e1_dp*t708*t246 - 0.12e2_dp*t77*t711*t682 &
    2045     18727793 :                                            + t745
    2046              :                   epsilon_c_unif1rhoa = e_c_u_01rhoa + t705*t82 + t708*t82 - &
    2047     18727793 :                                         t248 + t733*t80 + t736*t80 + t256
    2048     18727793 :                   t750 = 0.1e1_dp/t72
    2049     18727793 :                   t755 = 0.1e1_dp/t75
    2050              :                   phirhoarhoa = -t750*t682/0.9e1_dp + t257*chirhoarhoa/ &
    2051     18727793 :                                 0.3e1_dp - t755*t682/0.9e1_dp - t259*chirhoarhoa/0.3e1_dp
    2052     18727793 :                   phi1rhoa = phirhoa
    2053     18727793 :                   t763 = t90**2
    2054     18727793 :                   k_frhoarhoa = -0.2e1_dp/0.9e1_dp*t3/t262/t91*t763
    2055     18727793 :                   t767 = 0.1e1_dp/t94/t93
    2056     18727793 :                   t768 = k_frhoa**2
    2057     18727793 :                   k_s1rhoa = k_srhoa
    2058     18727793 :                   t775 = my_norm_drho*t105*t97
    2059     18727793 :                   t776 = t2*phirhoa
    2060     18727793 :                   t779 = t269*t273
    2061     18727793 :                   t785 = t269*t277*phirhoa/0.2e1_dp
    2062     18727793 :                   t789 = t2*k_srhoa
    2063     18727793 :                   t795 = t96/t272/k_s
    2064     18727793 :                   t798 = t273*t163
    2065     18727793 :                   t801 = t96*t798*k_srhoa/0.2e1_dp
    2066     18727793 :                   t812 = t96*t97*t518
    2067              :                   trhoarhoa = t775*t776*phi1rhoa + t779*t776*k_s1rhoa/ &
    2068              :                               0.2e1_dp + t785 - t269*t98*phirhoarhoa/0.2e1_dp + t779*t789 &
    2069              :                               *phi1rhoa/0.2e1_dp + t795*t789*k_s1rhoa + t801 - t96*t274* &
    2070              :                               (-t767*t768*t523/0.2e1_dp + t266*k_frhoarhoa*t7)/ &
    2071              :                               0.2e1_dp + t269*t277*phi1rhoa/0.2e1_dp + t96*t798*k_s1rhoa &
    2072     18727793 :                               /0.2e1_dp + t812
    2073              :                   t1rhoa = -t269*t98*phi1rhoa/0.2e1_dp - t96*t274*k_s1rhoa &
    2074     18727793 :                            /0.2e1_dp - t278/0.2e1_dp
    2075     18727793 :                   t820 = t101/t280/t108
    2076     18727793 :                   t821 = t107**2
    2077     18727793 :                   t822 = t289*t821
    2078     18727793 :                   t823 = epsilon_c_unif1rhoa*t100
    2079     18727793 :                   t825 = t285*phi1rhoa
    2080     18727793 :                   t828 = -t823*t105 + 0.3e1_dp*t102*t825
    2081     18727793 :                   t839 = 0.1e1_dp/t284/phi
    2082     18727793 :                   t840 = t839*phirhoa
    2083     18727793 :                   t851 = t101*t281
    2084              :                   Arhoarhoa = 0.2e1_dp*t820*t822*t828 - t101*t281*( &
    2085              :                               -epsilon_c_unifrhoarhoa*t100*t105 + 0.3e1_dp*t282*t825 + &
    2086              :                               0.3e1_dp*t823*t286 - 0.12e2_dp*t102*t840*phi1rhoa + &
    2087              :                               0.3e1_dp*t102*t285*phirhoarhoa)*t107 - t851*t289*t828* &
    2088     18727793 :                               t107
    2089     18727793 :                   A1rhoa = -t101*t281*t828*t107
    2090     18727793 :                   t858 = gamma_var*phi
    2091     18727793 :                   t865 = A1rhoa*t111
    2092     18727793 :                   t867 = 0.2e1_dp*t303*t1rhoa
    2093     18727793 :                   t868 = t865 + t867
    2094     18727793 :                   t876 = t865 + t867 + 0.2e1_dp*t314*A1rhoa + 0.4e1_dp*t318*t1rhoa
    2095              :                   t879 = 0.2e1_dp*t297*t298*t1rhoa + t101*t111*t868*t119 &
    2096     18727793 :                          - t310*t313*t876
    2097     18727793 :                   t880 = t879*t325
    2098     18727793 :                   t904 = t306*t119
    2099     18727793 :                   t908 = Arhoarhoa*t111
    2100     18727793 :                   t909 = Arhoa*t
    2101     18727793 :                   t911 = 0.2e1_dp*t909*t1rhoa
    2102     18727793 :                   t914 = 0.2e1_dp*A1rhoa*t*trhoa
    2103     18727793 :                   t917 = 0.2e1_dp*A*t1rhoa*trhoa
    2104     18727793 :                   t919 = 0.2e1_dp*t303*trhoarhoa
    2105     18727793 :                   t924 = t306*t312
    2106     18727793 :                   t936 = t113/t311/t118
    2107     18727793 :                   t944 = A*t317
    2108     18727793 :                   t953 = t115*t111
    2109              :                   t959 = t908 + t911 + t914 + t917 + t919 + 0.2e1_dp*A1rhoa*t116 &
    2110              :                          *Arhoa + 0.8e1_dp*t944*Arhoa*t1rhoa + 0.2e1_dp*t314* &
    2111              :                          Arhoarhoa + 0.8e1_dp*t944*trhoa*A1rhoa + 0.12e2_dp*t953* &
    2112     18727793 :                          trhoa*t1rhoa + 0.4e1_dp*t318*trhoarhoa
    2113              :                   t962 = 0.2e1_dp*t101*t1rhoa*t299 + 0.2e1_dp*t297*t868* &
    2114              :                          t119*trhoa - 0.2e1_dp*t297*t313*trhoa*t876 + 0.2e1_dp* &
    2115              :                          t297*t298*trhoarhoa + 0.2e1_dp*t297*t904*t1rhoa + t101* &
    2116              :                          t111*(t908 + t911 + t914 + t917 + t919)*t119 - t310*t924* &
    2117              :                          t876 - 0.2e1_dp*t297*t313*t321*t1rhoa - t310*t868*t312* &
    2118     18727793 :                          t321 + 0.2e1_dp*t310*t936*t321*t876 - t310*t313*t959
    2119     18727793 :                   t965 = t122**2
    2120     18727793 :                   t966 = 0.1e1_dp/t965
    2121     18727793 :                   t967 = t324*t966
    2122     18727793 :                   kf_arhoarhoa = -0.2e1_dp/0.9e1_dp*t124/t329/t125*t763
    2123     18727793 :                   ex_unif_a1rhoa = ex_unif_arhoa
    2124     18727793 :                   t985 = kf_arhoa**2
    2125     18727793 :                   s_a1rhoa = s_arhoa
    2126     18727793 :                   t998 = mu**2
    2127     18727793 :                   t999 = 0.1e1_dp/t344/t137*t998
    2128     18727793 :                   t1000 = t999*t133
    2129     18727793 :                   t1001 = s_arhoa*t135
    2130     18727793 :                   Fx_a1rhoa = 0.2e1_dp*t346*s_a*s_a1rhoa
    2131              : 
    2132              :                   e_ra_ra(ii) = e_ra_ra(ii) + &
    2133              :                                 scale_ex*(0.2e1_dp*ex_unif_a1rhoa*Fx_a + &
    2134              :                                           0.2e1_dp*ex_unif_a*Fx_a1rhoa + 0.2e1_dp*ex_unif_arhoa*Fx_a - &
    2135              :                                           0.3e1_dp/0.2e1_dp*my_rhoa*t7*kf_arhoarhoa*Fx_a + 0.2e1_dp* &
    2136              :                                           t350*Fx_a1rhoa + 0.2e1_dp*ex_unif_a*Fx_arhoa + 0.2e1_dp*my_rhoa &
    2137              :                                           *ex_unif_a1rhoa*Fx_arhoa + 0.2e1_dp*t140*(-0.8e1_dp*t1000 &
    2138              :                                                                        *t1001*s_a1rhoa + 0.2e1_dp*t346*s_a1rhoa*s_arhoa + 0.2e1_dp &
    2139              :                                                                               *t346*s_a*(my_norm_drhoa/t335/kf_a*t131*t985 + t337* &
    2140              :                                                                            t341*kf_arhoa - t337*t131*kf_arhoarhoa/0.2e1_dp + t130/ &
    2141              :                                                                         t340/my_rhoa)))/0.2e1_dp + scale_ec*(epsilon_c_unif1rhoa + &
    2142              :                                                                       0.3e1_dp*t293*t123*phi1rhoa + t110*t880 + epsilon_cGGArhoa + &
    2143              :                                                                     my_rho*(epsilon_c_unifrhoarhoa + 0.6e1_dp*t858*t294*phi1rhoa + &
    2144              :                                                                                   0.3e1_dp*t293*t880*phirhoa + 0.3e1_dp*t293*t123* &
    2145              :                                                                         phirhoarhoa + 0.3e1_dp*t293*t326*phi1rhoa + t110*t962*t325 &
    2146     18727793 :                                                                                                                   - t110*t967*t879))
    2147              : 
    2148     18727793 :                   chirhoarhob = 0.2e1_dp*t519
    2149     18727793 :                   rsrhoarhob = rsrhoarhoa
    2150     18727793 :                   t1031 = t176*t368*t191
    2151     18727793 :                   t1033 = alpha_1_1*rsrhob
    2152     18727793 :                   t1038 = rsrhoa*rsrhob
    2153     18727793 :                   t1050 = rsrhob*t564*rsrhoa
    2154              :                   e_c_u_0rhoarhob = -0.2e1_dp*t171*rsrhoarhob*t28 + t536* &
    2155              :                                     t1031 + t1033*t538 - 0.2e1_dp*t543*t192*t368 + t177*(-t549 &
    2156              :                                                                             *t1038/0.4e1_dp + t179*rsrhoarhob/0.2e1_dp + beta_2_1* &
    2157              :                                                                              rsrhoarhob + 0.3e1_dp/0.4e1_dp*t556*t1038 + 0.3e1_dp/ &
    2158              :                                                                               0.2e1_dp*t183*rsrhoarhob + t22*t561*t1050 + t22*t20* &
    2159              :                                                                            rsrhoarhob*t187 - t22*t20*t1050)*t191 + t578*t190*t580* &
    2160     18727793 :                                     t14*t368/0.2e1_dp
    2161     18727793 :                   t1069 = t199*t382*t212
    2162     18727793 :                   t1071 = alpha_1_2*rsrhob
    2163     18727793 :                   t1104 = t220*t396*t233
    2164     18727793 :                   t1106 = alpha_1_3*rsrhob
    2165              :                   frhoarhob = (0.4e1_dp/0.9e1_dp*t681*chirhoa*chirhob + &
    2166              :                                0.4e1_dp/0.3e1_dp*t71*chirhoarhob + 0.4e1_dp/0.9e1_dp*t687 &
    2167              :                                *chirhoa*chirhob - 0.4e1_dp/0.3e1_dp*t74*chirhoarhob)* &
    2168     18727793 :                               t69
    2169     18727793 :                   t1164 = t79*chirhoa*chirhob
    2170              :                   t1193 = -0.4e1_dp*t77*t245*chirhoarhob + (-0.2e1_dp*t194* &
    2171              :                                                             rsrhoarhob*t46 + t588*t1069 + t1071*t590 - 0.2e1_dp*t595* &
    2172              :                                                             t213*t382 + t200*(-t600*t1038/0.4e1_dp + t201*rsrhoarhob/ &
    2173              :                                                                           0.2e1_dp + beta_2_2*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t606* &
    2174              :                                                                         t1038 + 0.3e1_dp/0.2e1_dp*t205*rsrhoarhob + t40*t611*t1050 &
    2175              :                                                                             + t40*t38*rsrhoarhob*t187 - t40*t38*t1050)*t212 + t626 &
    2176              :                                                             *t211*t628*t34*t382/0.2e1_dp - e_c_u_0rhoarhob)*f*t80 + &
    2177              :                           t249*frhob*t80 + 0.4e1_dp*t250*t415 + t410*frhoa*t80 + &
    2178              :                           t84*frhoarhob*t80 + 0.4e1_dp*t252*t415 + 0.4e1_dp*t411* &
    2179              :                           t254 + 0.4e1_dp*t413*t254 + 0.12e2_dp*t85*t1164 + 0.4e1_dp* &
    2180     18727793 :                           t85*t244*chirhoarhob
    2181              :                   epsilon_c_unifrhoarhob = e_c_u_0rhoarhob + (0.2e1_dp*t215* &
    2182              :                                                               rsrhoarhob*t64 - t636*t1104 - t1106*t638 + 0.2e1_dp*t643* &
    2183              :                                                               t234*t396 - t221*(-t648*t1038/0.4e1_dp + t222*rsrhoarhob/ &
    2184              :                                                                           0.2e1_dp + beta_2_3*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t654* &
    2185              :                                                                         t1038 + 0.3e1_dp/0.2e1_dp*t226*rsrhoarhob + t58*t659*t1050 &
    2186              :                                                                             + t58*t56*rsrhoarhob*t187 - t58*t56*t1050)*t233 - t674 &
    2187              :                                                               *t232*t676*t52*t396/0.2e1_dp)*f*t82 + alpha_crhoa* &
    2188              :                                            frhob*t82 - 0.4e1_dp*t240*t407 + alpha_crhob*frhoa*t82 + &
    2189              :                                            alpha_c*frhoarhob*t82 - 0.4e1_dp*t242*t407 - 0.4e1_dp*t403 &
    2190              :                                            *t246 - 0.4e1_dp*t405*t246 - 0.12e2_dp*t77*t78*t1164 + &
    2191     18727793 :                                            t1193
    2192              :                   phirhoarhob = -t750*chirhoa*chirhob/0.9e1_dp + t257* &
    2193              :                                 chirhoarhob/0.3e1_dp - t755*chirhoa*chirhob/0.9e1_dp - t259 &
    2194     18727793 :                                 *chirhoarhob/0.3e1_dp
    2195     18727793 :                   k_frhoarhob = k_frhoarhoa
    2196     18727793 :                   t1228 = t269*t277*phirhob/0.2e1_dp
    2197     18727793 :                   t1231 = t96*t798*k_srhob/0.2e1_dp
    2198              :                   trhoarhob = t775*t776*phirhob + t779*t776*k_srhob/ &
    2199              :                               0.2e1_dp + t785 - t269*t98*phirhoarhob/0.2e1_dp + t779*t789 &
    2200              :                               *phirhob/0.2e1_dp + t795*t789*k_srhob + t801 - t96*t274*( &
    2201              :                               -t767*k_frhoa*t523*k_frhob/0.2e1_dp + t266*k_frhoarhob* &
    2202     18727793 :                               t7)/0.2e1_dp + t1228 + t1231 + t812
    2203              :                   Arhoarhob = 0.2e1_dp*t820*t822*t432 - t101*t281*( &
    2204              :                               -epsilon_c_unifrhoarhob*t100*t105 + 0.3e1_dp*t282*t429 + &
    2205              :                               0.3e1_dp*t427*t286 - 0.12e2_dp*t102*t840*phirhob + &
    2206              :                               0.3e1_dp*t102*t285*phirhoarhob)*t107 - t851*t289*t432* &
    2207     18727793 :                               t107
    2208     18727793 :                   t1269 = t445*t119
    2209     18727793 :                   t1283 = Arhoarhob*t111
    2210     18727793 :                   t1285 = 0.2e1_dp*t909*trhob
    2211     18727793 :                   t1286 = Arhob*t
    2212     18727793 :                   t1288 = 0.2e1_dp*t1286*trhoa
    2213     18727793 :                   t1291 = 0.2e1_dp*A*trhob*trhoa
    2214     18727793 :                   t1293 = 0.2e1_dp*t303*trhoarhob
    2215     18727793 :                   t1304 = t445*t312
    2216              :                   t1327 = t1283 + t1285 + t1288 + t1291 + t1293 + 0.2e1_dp*Arhob* &
    2217              :                           t116*Arhoa + 0.8e1_dp*t944*Arhoa*trhob + 0.2e1_dp*t314* &
    2218              :                           Arhoarhob + 0.8e1_dp*t944*trhoa*Arhob + 0.12e2_dp*t953* &
    2219     18727793 :                           trhoa*trhob + 0.4e1_dp*t318*trhoarhob
    2220              :                   t1330 = 0.2e1_dp*t101*trhob*t299 + 0.2e1_dp*t297*t1269* &
    2221              :                           trhoa - 0.2e1_dp*t297*t313*trhoa*t453 + 0.2e1_dp*t297* &
    2222              :                           t298*trhoarhob + 0.2e1_dp*t297*t904*trhob + t101*t111*( &
    2223              :                           t1283 + t1285 + t1288 + t1291 + t1293)*t119 - t310*t924*t453 - &
    2224              :                           0.2e1_dp*t297*t313*t321*trhob - t310*t1304*t321 + &
    2225     18727793 :                           0.2e1_dp*t310*t936*t321*t453 - t310*t313*t1327
    2226              : 
    2227              :                   e_ra_rb(ii) = e_ra_rb(ii) + &
    2228              :                                 scale_ec*(epsilon_cGGArhob + epsilon_cGGArhoa + &
    2229              :                                           my_rho*(epsilon_c_unifrhoarhob + 0.6e1_dp*t858*t294*phirhob + &
    2230              :                                                   0.3e1_dp*t293*t457*phirhoa + 0.3e1_dp*t293*t123* &
    2231              :                                                   phirhoarhob + 0.3e1_dp*t293*t326*phirhob + t110*t1330*t325 &
    2232     18727793 :                                                   - t110*t967*t456))
    2233              : 
    2234     18727793 :                   chirhobrhob = 0.2e1_dp*t163 + 0.2e1_dp*t519
    2235     18727793 :                   rsrhobrhob = rsrhoarhob
    2236     18727793 :                   t1342 = t368**2
    2237     18727793 :                   t1346 = rsrhob**2
    2238              :                   e_c_u_0rhobrhob = -0.2e1_dp*t171*rsrhobrhob*t28 + 0.2e1_dp* &
    2239              :                                     t1033*t1031 - 0.2e1_dp*t543*t1342*t191 + t177*(-t549* &
    2240              :                                                                              t1346/0.4e1_dp + t179*rsrhobrhob/0.2e1_dp + beta_2_1* &
    2241              :                                                                              rsrhobrhob + 0.3e1_dp/0.4e1_dp*t556*t1346 + 0.3e1_dp/ &
    2242              :                                                                           0.2e1_dp*t183*rsrhobrhob + t22*t561*t1346*t564 + t22*t20 &
    2243              :                                                                                *rsrhobrhob*t187 - t22*t20*t1346*t564)*t191 + t578* &
    2244     18727793 :                                     t1342*t580*t14/0.2e1_dp
    2245     18727793 :                   e_c_u_01rhob = e_c_u_0rhob
    2246     18727793 :                   t1377 = t382**2
    2247     18727793 :                   t1411 = t396**2
    2248     18727793 :                   alpha_c1rhob = alpha_crhob
    2249     18727793 :                   t1440 = chirhob**2
    2250              :                   frhobrhob = (0.4e1_dp/0.9e1_dp*t681*t1440 + 0.4e1_dp/ &
    2251              :                                0.3e1_dp*t71*chirhobrhob + 0.4e1_dp/0.9e1_dp*t687*t1440 - &
    2252     18727793 :                                0.4e1_dp/0.3e1_dp*t74*chirhobrhob)*t69
    2253     18727793 :                   f1rhob = frhob
    2254     18727793 :                   t1462 = alpha_c1rhob*f
    2255     18727793 :                   t1465 = alpha_c*f1rhob
    2256     18727793 :                   t1482 = e_c_u_1rhob - e_c_u_01rhob
    2257     18727793 :                   t1489 = t1482*f
    2258     18727793 :                   t1492 = t84*f1rhob
    2259              :                   t1501 = -0.4e1_dp*t77*t245*chirhobrhob + (-0.2e1_dp*t194* &
    2260              :                                                             rsrhobrhob*t46 + 0.2e1_dp*t1071*t1069 - 0.2e1_dp*t595* &
    2261              :                                                             t1377*t212 + t200*(-t600*t1346/0.4e1_dp + t201*rsrhobrhob &
    2262              :                                                                          /0.2e1_dp + beta_2_2*rsrhobrhob + 0.3e1_dp/0.4e1_dp*t606* &
    2263              :                                                                         t1346 + 0.3e1_dp/0.2e1_dp*t205*rsrhobrhob + t40*t611*t1346 &
    2264              :                                                                              *t564 + t40*t38*rsrhobrhob*t187 - t40*t38*t1346*t564) &
    2265              :                                                             *t212 + t626*t1377*t628*t34/0.2e1_dp - e_c_u_0rhobrhob)*f &
    2266              :                           *t80 + t410*f1rhob*t80 + 0.4e1_dp*t411*t415 + t1482* &
    2267              :                           frhob*t80 + t84*frhobrhob*t80 + 0.4e1_dp*t413*t415 + &
    2268              :                           0.4e1_dp*t1489*t415 + 0.4e1_dp*t1492*t415 + 0.12e2_dp*t85 &
    2269     18727793 :                           *t79*t1440 + 0.4e1_dp*t85*t244*chirhobrhob
    2270              :                   epsilon_c_unifrhobrhob = e_c_u_0rhobrhob + (0.2e1_dp*t215* &
    2271              :                                                               rsrhobrhob*t64 - 0.2e1_dp*t1106*t1104 + 0.2e1_dp*t643* &
    2272              :                                                               t1411*t233 - t221*(-t648*t1346/0.4e1_dp + t222*rsrhobrhob &
    2273              :                                                                          /0.2e1_dp + beta_2_3*rsrhobrhob + 0.3e1_dp/0.4e1_dp*t654* &
    2274              :                                                                         t1346 + 0.3e1_dp/0.2e1_dp*t226*rsrhobrhob + t58*t659*t1346 &
    2275              :                                                                              *t564 + t58*t56*rsrhobrhob*t187 - t58*t56*t1346*t564) &
    2276              :                                                               *t233 - t674*t1411*t676*t52/0.2e1_dp)*f*t82 + &
    2277              :                                            alpha_crhob*f1rhob*t82 - 0.4e1_dp*t403*t407 + alpha_c1rhob* &
    2278              :                                            frhob*t82 + alpha_c*frhobrhob*t82 - 0.4e1_dp*t405*t407 - &
    2279              :                                            0.4e1_dp*t1462*t407 - 0.4e1_dp*t1465*t407 - 0.12e2_dp*t77 &
    2280     18727793 :                                            *t711*t1440 + t1501
    2281              :                   epsilon_c_unif1rhob = e_c_u_01rhob + t1462*t82 + t1465*t82 - &
    2282     18727793 :                                         t409 + t1489*t80 + t1492*t80 + t417
    2283              :                   phirhobrhob = -t750*t1440/0.9e1_dp + t257*chirhobrhob/ &
    2284     18727793 :                                 0.3e1_dp - t755*t1440/0.9e1_dp - t259*chirhobrhob/0.3e1_dp
    2285     18727793 :                   phi1rhob = phirhob
    2286     18727793 :                   t1514 = k_frhob**2
    2287     18727793 :                   k_s1rhob = k_srhob
    2288     18727793 :                   t1520 = t2*phirhob
    2289     18727793 :                   t1529 = t2*k_srhob
    2290              :                   trhobrhob = t775*t1520*phi1rhob + t779*t1520*k_s1rhob/ &
    2291              :                               0.2e1_dp + t1228 - t269*t98*phirhobrhob/0.2e1_dp + t779* &
    2292              :                               t1529*phi1rhob/0.2e1_dp + t795*t1529*k_s1rhob + t1231 - t96 &
    2293              :                               *t274*(-t767*t1514*t523/0.2e1_dp + t266*k_frhoarhob*t7) &
    2294              :                               /0.2e1_dp + t269*t277*phi1rhob/0.2e1_dp + t96*t798* &
    2295     18727793 :                               k_s1rhob/0.2e1_dp + t812
    2296              :                   t1rhob = -t269*t98*phi1rhob/0.2e1_dp - t96*t274*k_s1rhob &
    2297     18727793 :                            /0.2e1_dp - t278/0.2e1_dp
    2298     18727793 :                   t1550 = epsilon_c_unif1rhob*t100
    2299     18727793 :                   t1552 = t285*phi1rhob
    2300     18727793 :                   t1555 = -t1550*t105 + 0.3e1_dp*t102*t1552
    2301              :                   Arhobrhob = 0.2e1_dp*t820*t432*t821*t1555 - t101*t281* &
    2302              :                               (-epsilon_c_unifrhobrhob*t100*t105 + 0.3e1_dp*t427*t1552 + &
    2303              :                                0.3e1_dp*t1550*t429 - 0.12e2_dp*t102*t839*phirhob* &
    2304              :                                phi1rhob + 0.3e1_dp*t102*t285*phirhobrhob)*t107 - t851* &
    2305     18727793 :                               t432*t1555*t107
    2306     18727793 :                   A1rhob = -t101*t281*t1555*t107
    2307     18727793 :                   t1588 = A1rhob*t111
    2308     18727793 :                   t1590 = 0.2e1_dp*t303*t1rhob
    2309     18727793 :                   t1591 = t1588 + t1590
    2310              :                   t1599 = t1588 + t1590 + 0.2e1_dp*t314*A1rhob + 0.4e1_dp*t318 &
    2311     18727793 :                           *t1rhob
    2312              :                   t1602 = 0.2e1_dp*t297*t298*t1rhob + t101*t111*t1591* &
    2313     18727793 :                           t119 - t310*t313*t1599
    2314     18727793 :                   t1603 = t1602*t325
    2315     18727793 :                   t1630 = Arhobrhob*t111
    2316     18727793 :                   t1632 = 0.2e1_dp*t1286*t1rhob
    2317     18727793 :                   t1635 = 0.2e1_dp*A1rhob*t*trhob
    2318     18727793 :                   t1638 = 0.2e1_dp*A*t1rhob*trhob
    2319     18727793 :                   t1640 = 0.2e1_dp*t303*trhobrhob
    2320              :                   t1674 = t1630 + t1632 + t1635 + t1638 + t1640 + 0.2e1_dp*A1rhob &
    2321              :                           *t116*Arhob + 0.8e1_dp*t944*Arhob*t1rhob + 0.2e1_dp*t314 &
    2322              :                           *Arhobrhob + 0.8e1_dp*t944*trhob*A1rhob + 0.12e2_dp*t953* &
    2323     18727793 :                           trhob*t1rhob + 0.4e1_dp*t318*trhobrhob
    2324              :                   t1677 = 0.2e1_dp*t101*t1rhob*t439 + 0.2e1_dp*t297*t1591 &
    2325              :                           *t119*trhob - 0.2e1_dp*t297*t313*trhob*t1599 + 0.2e1_dp* &
    2326              :                           t297*t298*trhobrhob + 0.2e1_dp*t297*t1269*t1rhob + t101* &
    2327              :                           t111*(t1630 + t1632 + t1635 + t1638 + t1640)*t119 - t310* &
    2328              :                           t1304*t1599 - 0.2e1_dp*t297*t313*t453*t1rhob - t310* &
    2329              :                           t1591*t312*t453 + 0.2e1_dp*t310*t936*t453*t1599 - t310* &
    2330     18727793 :                           t313*t1674
    2331     18727793 :                   t1680 = t456*t966
    2332     18727793 :                   kf_brhobrhob = -0.2e1_dp/0.9e1_dp*t124/t460/t142*t763
    2333     18727793 :                   ex_unif_b1rhob = ex_unif_brhob
    2334     18727793 :                   t1698 = kf_brhob**2
    2335     18727793 :                   s_b1rhob = s_brhob
    2336     18727793 :                   t1711 = 0.1e1_dp/t475/t153*t998
    2337     18727793 :                   t1712 = t1711*t150
    2338     18727793 :                   t1713 = s_brhob*t135
    2339     18727793 :                   Fx_b1rhob = 0.2e1_dp*t477*s_b*s_b1rhob
    2340              : 
    2341              :                   e_rb_rb(ii) = e_rb_rb(ii) + &
    2342              :                                 scale_ex*(0.2e1_dp*ex_unif_b1rhob*Fx_b + &
    2343              :                                           0.2e1_dp*ex_unif_b*Fx_b1rhob + 0.2e1_dp*ex_unif_brhob*Fx_b - &
    2344              :                                           0.3e1_dp/0.2e1_dp*my_rhob*t7*kf_brhobrhob*Fx_b + 0.2e1_dp* &
    2345              :                                           t481*Fx_b1rhob + 0.2e1_dp*ex_unif_b*Fx_brhob + 0.2e1_dp*my_rhob &
    2346              :                                           *ex_unif_b1rhob*Fx_brhob + 0.2e1_dp*t156*(-0.8e1_dp*t1712 &
    2347              :                                                                        *t1713*s_b1rhob + 0.2e1_dp*t477*s_b1rhob*s_brhob + 0.2e1_dp &
    2348              :                                                                              *t477*s_b*(my_norm_drhob/t466/kf_b*t148*t1698 + t468* &
    2349              :                                                                            t472*kf_brhob - t468*t148*kf_brhobrhob/0.2e1_dp + t147/ &
    2350              :                                                                         t471/my_rhob)))/0.2e1_dp + scale_ec*(epsilon_c_unif1rhob + &
    2351              :                                                                        0.3e1_dp*t293*t123*phi1rhob + t110*t1603 + epsilon_cGGArhob &
    2352              :                                                                     + my_rho*(epsilon_c_unifrhobrhob + 0.6e1_dp*t858*t436*phi1rhob &
    2353              :                                                                                + 0.3e1_dp*t293*t1603*phirhob + 0.3e1_dp*t293*t123* &
    2354              :                                                                            phirhobrhob + 0.3e1_dp*t293*t457*phi1rhob + t110*t1677* &
    2355     18727793 :                                                                                                            t325 - t110*t1680*t1602))
    2356     18727793 :                   t1739 = t268*t97
    2357     18727793 :                   t1741 = t95*t273
    2358     18727793 :                   t1743 = t488*t163
    2359              :                   trhoanorm_drho = -t1739*t776/0.2e1_dp - t1741*t789/ &
    2360     18727793 :                                    0.2e1_dp - t1743/0.2e1_dp
    2361     18727793 :                   t1748 = t101*tnorm_drho
    2362     18727793 :                   t1765 = t909*tnorm_drho
    2363     18727793 :                   t1766 = t494*trhoa
    2364     18727793 :                   t1767 = t303*trhoanorm_drho
    2365              :                   t1801 = 0.2e1_dp*t1748*t299 + 0.4e1_dp*t310*t494*t119* &
    2366              :                           trhoa - 0.2e1_dp*t297*t313*trhoa*t502 + 0.2e1_dp*t297* &
    2367              :                           t298*trhoanorm_drho + 0.2e1_dp*t297*t904*tnorm_drho + t101* &
    2368              :                           t111*(0.2e1_dp*t1765 + 0.2e1_dp*t1766 + 0.2e1_dp*t1767)* &
    2369              :                           t119 - t310*t924*t502 - 0.2e1_dp*t297*t313*t321* &
    2370              :                           tnorm_drho - 0.2e1_dp*t493*t494*t312*t321 + 0.2e1_dp*t310 &
    2371              :                           *t936*t321*t502 - t310*t313*(0.2e1_dp*t1765 + 0.2e1_dp* &
    2372              :                                                        t1766 + 0.2e1_dp*t1767 + 0.8e1_dp*t944*Arhoa*tnorm_drho + &
    2373              :                                                        0.12e2_dp*t953*trhoa*tnorm_drho + 0.4e1_dp*t318* &
    2374     18727793 :                                                        trhoanorm_drho)
    2375              : 
    2376              :                   e_ra_ndr(ii) = e_ra_ndr(ii) + &
    2377              :                                  scale_ec*(Hnorm_drho + my_rho*(0.3e1_dp* &
    2378     18727793 :                                                                 t293*t506*phirhoa + t110*t1801*t325 - t110*t967*t505))
    2379              : 
    2380              :                   trhobnorm_drho = -t1739*t1520/0.2e1_dp - t1741*t1529/ &
    2381     18727793 :                                    0.2e1_dp - t1743/0.2e1_dp
    2382     18727793 :                   t1829 = t1286*tnorm_drho
    2383     18727793 :                   t1830 = t494*trhob
    2384     18727793 :                   t1831 = t303*trhobnorm_drho
    2385              :                   t1865 = 0.2e1_dp*t1748*t439 + 0.4e1_dp*t310*t494*t119* &
    2386              :                           trhob - 0.2e1_dp*t297*t313*trhob*t502 + 0.2e1_dp*t297* &
    2387              :                           t298*trhobnorm_drho + 0.2e1_dp*t297*t1269*tnorm_drho + t101 &
    2388              :                           *t111*(0.2e1_dp*t1829 + 0.2e1_dp*t1830 + 0.2e1_dp*t1831)* &
    2389              :                           t119 - t310*t1304*t502 - 0.2e1_dp*t297*t313*t453* &
    2390              :                           tnorm_drho - 0.2e1_dp*t493*t494*t312*t453 + 0.2e1_dp*t310 &
    2391              :                           *t936*t453*t502 - t310*t313*(0.2e1_dp*t1829 + 0.2e1_dp* &
    2392              :                                                        t1830 + 0.2e1_dp*t1831 + 0.8e1_dp*t944*Arhob*tnorm_drho + &
    2393              :                                                        0.12e2_dp*t953*trhob*tnorm_drho + 0.4e1_dp*t318* &
    2394     18727793 :                                                        trhobnorm_drho)
    2395              : 
    2396              :                   e_rb_ndr(ii) = e_rb_ndr(ii) + &
    2397              :                                  scale_ec*(Hnorm_drho + my_rho*(0.3e1_dp* &
    2398     18727793 :                                                                 t293*t506*phirhob + t110*t1865*t325 - t110*t1680*t505))
    2399              : 
    2400     18727793 :                   t1871 = tnorm_drho**2
    2401     18727793 :                   t1876 = A*t1871
    2402     18727793 :                   t1888 = t502**2
    2403     18727793 :                   t1901 = t505**2
    2404              : 
    2405              :                   e_ndr_ndr(ii) = e_ndr_ndr(ii) + &
    2406              :                                   scale_ec*my_rho*(t110*(0.2e1_dp* &
    2407              :                                                          t101*t1871*t113*t119 + 0.10e2_dp*t310*t1876*t119 - &
    2408              :                                                          0.4e1_dp*t297*t313*tnorm_drho*t502 - 0.4e1_dp*t493*t494 &
    2409              :                                                          *t312*t502 + 0.2e1_dp*t310*t936*t1888 - t310*t313*( &
    2410              :                                                          0.2e1_dp*t1876 + 0.12e2_dp*t953*t1871))*t325 - t110*t1901 &
    2411     18727793 :                                                    *t966)
    2412              : 
    2413              :                   e_ra_ndra(ii) = e_ra_ndra(ii) + &
    2414              :                                   scale_ex*(0.2e1_dp*ex_unif_a* &
    2415              :                                             Fx_anorm_drhoa + 0.2e1_dp*t350*Fx_anorm_drhoa + 0.2e1_dp*t140 &
    2416              :                                             *(-0.8e1_dp*t1000*t1001*s_anorm_drhoa + 0.2e1_dp*t346* &
    2417              :                                               s_anorm_drhoa*s_arhoa + 0.2e1_dp*t346*s_a*(-t336*t131* &
    2418     18727793 :                                                                                   kf_arhoa/0.2e1_dp - t129*t341/0.2e1_dp)))/0.2e1_dp
    2419              : 
    2420     18727793 :                   t1922 = s_anorm_drhoa**2
    2421              : 
    2422              :                   e_ndra_ndra(ii) = e_ndra_ndra(ii) + &
    2423              :                                     scale_ex*t140*(-0.8e1_dp*t999* &
    2424     18727793 :                                                    t133*t1922*t135 + 0.2e1_dp*t346*t1922)
    2425              : 
    2426              :                   e_rb_ndrb(ii) = e_rb_ndrb(ii) + &
    2427              :                                   scale_ex*(0.2e1_dp*ex_unif_b* &
    2428              :                                             Fx_bnorm_drhob + 0.2e1_dp*t481*Fx_bnorm_drhob + 0.2e1_dp*t156 &
    2429              :                                             *(-0.8e1_dp*t1712*t1713*s_bnorm_drhob + 0.2e1_dp*t477* &
    2430              :                                               s_bnorm_drhob*s_brhob + 0.2e1_dp*t477*s_b*(-t467*t148* &
    2431     18727793 :                                                                                   kf_brhob/0.2e1_dp - t146*t472/0.2e1_dp)))/0.2e1_dp
    2432              : 
    2433     18727793 :                   t1949 = s_bnorm_drhob**2
    2434              :                   e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii) + &
    2435              :                                     scale_ex*t156*(-0.8e1_dp*t1711* &
    2436     18727793 :                                                    t150*t1949*t135 + 0.2e1_dp*t477*t1949)
    2437              :                END IF
    2438              :             END IF
    2439              :          END DO
    2440              : !$OMP END DO
    2441              :       END SELECT
    2442        19263 :    END SUBROUTINE pbe_lsd_calc
    2443              : 
    2444              : END MODULE xc_pbe
        

Generated by: LCOV version 2.0-1