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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief calculates the 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      242163 :    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       80721 :       CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
      79       80721 :       CALL section_vals_val_get(pbe_params, "scale_x", r_val=sx)
      80       80721 :       CALL section_vals_val_get(pbe_params, "scale_c", r_val=sc)
      81             : 
      82       80368 :       SELECT CASE (param)
      83             :       CASE (xc_pbe_orig)
      84       80368 :          CALL cite_reference(Perdew1996)
      85       80368 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
      86       52136 :             IF (PRESENT(reference)) THEN
      87             :                reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "// &
      88         367 :                            "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996) {spin unpolarized}"
      89             :             END IF
      90       52136 :             IF (PRESENT(shortform)) THEN
      91         367 :                shortform = "PBE Perdew-Burke-Ernzerhof xc functional (unpolarized)"
      92             :             END IF
      93             :          ELSE
      94       28232 :             IF (PRESENT(reference)) THEN
      95             :                WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") &
      96         167 :                   "J.P.Perdew, K.Burke, M.Ernzerhof, ", &
      97         167 :                   "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)", &
      98         334 :                   sx, sc, "{spin unpolarized}"
      99             :             END IF
     100       28232 :             IF (PRESENT(shortform)) THEN
     101             :                WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") &
     102         167 :                   "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       80721 :          CPABORT("Unsupported parametrization")
     154             :       END SELECT
     155       80721 :       IF (PRESENT(needs)) THEN
     156       80180 :          needs%rho = .TRUE.
     157       80180 :          needs%norm_drho = .TRUE.
     158             :       END IF
     159       80721 :       IF (PRESENT(max_deriv)) max_deriv = 3
     160             : 
     161       80721 :    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       55044 :    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       18348 :       CALL section_vals_val_get(pbe_params, "parametrization", i_val=param)
     183       18348 :       CALL section_vals_val_get(pbe_params, "scale_x", r_val=sx)
     184       18348 :       CALL section_vals_val_get(pbe_params, "scale_c", r_val=sc)
     185             : 
     186       18176 :       SELECT CASE (param)
     187             :       CASE (xc_pbe_orig)
     188       18176 :          CALL cite_reference(Perdew1996)
     189       18176 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
     190       13172 :             IF (PRESENT(reference)) THEN
     191             :                reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "// &
     192         245 :                            "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996) {spin polarized}"
     193             :             END IF
     194       13172 :             IF (PRESENT(shortform)) THEN
     195         245 :                shortform = "PBE xc functional (spin polarized)"
     196             :             END IF
     197             :          ELSE
     198        5004 :             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        5004 :             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       18348 :          CPABORT("Unsupported parametrization")
     258             :       END SELECT
     259       18348 :       IF (PRESENT(needs)) THEN
     260       18051 :          needs%rho_spin = .TRUE.
     261       18051 :          needs%norm_drho_spin = .TRUE.
     262       18051 :          needs%norm_drho = .TRUE.
     263             :       END IF
     264       18348 :       IF (PRESENT(max_deriv)) max_deriv = 2
     265             : 
     266       18348 :    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      421890 :    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       84378 :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
     291       84378 :          e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
     292       84378 :          e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
     293             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     294             : 
     295       84378 :       CALL timeset(routineN, handle)
     296             : 
     297             :       CALL xc_rho_set_get(rho_set, rho=rho, &
     298       84378 :                           norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
     299       84378 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     300             : 
     301       84378 :       dummy => rho
     302             : 
     303       84378 :       e_0 => dummy
     304       84378 :       e_rho => dummy
     305       84378 :       e_ndrho => dummy
     306       84378 :       e_rho_rho => dummy
     307       84378 :       e_ndrho_rho => dummy
     308       84378 :       e_ndrho_ndrho => dummy
     309       84378 :       e_rho_rho_rho => dummy
     310       84378 :       e_ndrho_rho_rho => dummy
     311       84378 :       e_ndrho_ndrho_rho => dummy
     312       84378 :       e_ndrho_ndrho_ndrho => dummy
     313             : 
     314       84378 :       IF (grad_deriv >= 0) THEN
     315             :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     316       84378 :                                          allocate_deriv=.TRUE.)
     317       84378 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     318             :       END IF
     319       84378 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     320             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     321       82388 :                                          allocate_deriv=.TRUE.)
     322       82388 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     323             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     324       82388 :                                          allocate_deriv=.TRUE.)
     325       82388 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     326             :       END IF
     327       84378 :       IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
     328             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     329       11706 :                                          allocate_deriv=.TRUE.)
     330       11706 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     331             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
     332       11706 :                                          allocate_deriv=.TRUE.)
     333       11706 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
     334             :          deriv => xc_dset_get_derivative(deriv_set, &
     335       11706 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     336       11706 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     337             :       END IF
     338       84378 :       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       84378 :       IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
     353           0 :          CPABORT("derivatives bigger than 3 not implemented")
     354             :       END IF
     355             : 
     356       84378 :       CALL section_vals_val_get(pbe_params, "scale_c", r_val=scale_ec)
     357       84378 :       CALL section_vals_val_get(pbe_params, "scale_x", r_val=scale_ex)
     358       84378 :       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       84378 : !$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       84378 :       CALL timestop(handle)
     377       84378 :    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       84378 :    SUBROUTINE pbe_lda_calc(rho, norm_drho, &
     403       84378 :                            e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
     404       84378 :                            e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
     405       84378 :                            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       84378 :       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       84378 :          CPABORT("Unsupported parametrization")
     487             :       END SELECT
     488             : 
     489       84378 :       gamma_var = (0.1e1_dp - LOG(0.2e1_dp))/pi**2
     490       84378 :       p_1 = 0.10e1_dp
     491       84378 :       A_1 = 0.31091e-1_dp
     492       84378 :       alpha_1_1 = 0.21370e0_dp
     493       84378 :       beta_1_1 = 0.75957e1_dp
     494       84378 :       beta_2_1 = 0.35876e1_dp
     495       84378 :       beta_3_1 = 0.16382e1_dp
     496       84378 :       beta_4_1 = 0.49294e0_dp
     497       84378 :       p_2 = 0.10e1_dp
     498       84378 :       t1 = 3**(0.1e1_dp/0.3e1_dp)
     499       84378 :       t2 = 4**(0.1e1_dp/0.3e1_dp)
     500       84378 :       t3 = t2**2
     501       84378 :       t4 = t1*t3
     502       84378 :       t5 = 0.1e1_dp/pi
     503       84378 :       t69 = pi**2
     504       84378 :       t627 = 0.1e1_dp/t69/pi
     505             : 
     506       84378 :       SELECT CASE (grad_deriv)
     507             :       CASE default
     508             : !$OMP DO
     509             :          DO ii = 1, npoints
     510  1799645588 :             my_rho = rho(ii)
     511  1799645588 :             IF (my_rho > epsilon_rho) THEN
     512  1669325983 :                my_norm_drho = norm_drho(ii)
     513             : 
     514  1669325983 :                t6 = 0.1e1_dp/my_rho
     515  1669325983 :                t7 = t5*t6
     516  1669325983 :                t8 = t7**(0.1e1_dp/0.3e1_dp)
     517  1669325983 :                rs = t4*t8/0.4e1_dp
     518  1669325983 :                t11 = 0.1e1_dp + alpha_1_1*rs
     519  1669325983 :                t13 = 0.1e1_dp/A_1
     520  1669325983 :                t14 = SQRT(rs)
     521  1669325983 :                t17 = t14*rs
     522  1669325983 :                t19 = p_1 + 0.1e1_dp
     523  1669325983 :                t20 = rs**t19
     524  1669325983 :                t21 = beta_4_1*t20
     525  1669325983 :                t22 = beta_1_1*t14 + beta_2_1*rs + beta_3_1*t17 + t21
     526  1669325983 :                t26 = 0.1e1_dp + t13/t22/0.2e1_dp
     527  1669325983 :                t27 = LOG(t26)
     528  1669325983 :                e_c_u_0 = -0.2e1_dp*A_1*t11*t27
     529  1669325983 :                t65 = 2**(0.1e1_dp/0.3e1_dp)
     530  1669325983 :                t70 = t69*my_rho
     531  1669325983 :                t71 = t70**(0.1e1_dp/0.3e1_dp)
     532  1669325983 :                k_f = t1*t71
     533  1669325983 :                t72 = k_f*t5
     534  1669325983 :                t73 = SQRT(t72)
     535  1669325983 :                k_s = 0.2e1_dp*t73
     536  1669325983 :                t74 = 0.1e1_dp/k_s
     537  1669325983 :                t75 = my_norm_drho*t74
     538  1669325983 :                t = t75*t6/0.2e1_dp
     539  1669325983 :                t77 = 0.1e1_dp/gamma_var
     540  1669325983 :                t78 = beta*t77
     541  1669325983 :                t80 = EXP(-e_c_u_0*t77)
     542  1669325983 :                t81 = -0.1e1_dp + t80
     543  1669325983 :                A = t78/t81
     544  1669325983 :                t83 = t**2
     545  1669325983 :                t84 = A*t83
     546  1669325983 :                t85 = 0.1e1_dp + t84
     547  1669325983 :                t86 = t83*t85
     548  1669325983 :                t87 = A**2
     549  1669325983 :                t88 = t83**2
     550  1669325983 :                t90 = 0.1e1_dp + t84 + t87*t88
     551  1669325983 :                t91 = 0.1e1_dp/t90
     552  1669325983 :                t94 = 0.1e1_dp + t78*t86*t91
     553  1669325983 :                t95 = LOG(t94)
     554  1669325983 :                epsilon_cGGA = e_c_u_0 + gamma_var*t95
     555  1669325983 :                kf = k_f
     556  1669325983 :                ex_unif = -0.3e1_dp/0.4e1_dp*t5*kf
     557  1669325983 :                t98 = 0.1e1_dp/kf
     558  1669325983 :                t99 = my_norm_drho*t98
     559  1669325983 :                s = t99*t6/0.2e1_dp
     560  1669325983 :                t101 = s**2
     561  1669325983 :                t103 = 0.1e1_dp/kappa
     562  1669325983 :                t105 = 0.1e1_dp + mu*t101*t103
     563  1669325983 :                Fx = 0.1e1_dp + kappa - kappa/t105
     564  1669325983 :                t108 = my_rho*ex_unif
     565             : 
     566  1669325983 :                IF (grad_deriv >= 0) THEN
     567             :                   e_0(ii) = e_0(ii) + &
     568  1669325983 :                             scale_ex*t108*Fx + scale_ec*my_rho*epsilon_cGGA
     569             :                END IF
     570             : 
     571  1669325983 :                t111 = t8**2
     572  1669325983 :                t113 = 0.1e1_dp/t111*t5
     573  1669325983 :                t114 = my_rho**2
     574  1669325983 :                t115 = 0.1e1_dp/t114
     575  1669325983 :                rsrho = -t4*t113*t115/0.12e2_dp
     576  1669325983 :                t119 = A_1*alpha_1_1
     577  1669325983 :                t123 = t22**2
     578  1669325983 :                t124 = 0.1e1_dp/t123
     579  1669325983 :                t125 = t11*t124
     580  1669325983 :                t126 = 0.1e1_dp/t14
     581  1669325983 :                t127 = beta_1_1*t126
     582  1669325983 :                t131 = beta_3_1*t14
     583  1669325983 :                t135 = 0.1e1_dp/rs
     584             :                t138 = t127*rsrho/0.2e1_dp + beta_2_1*rsrho + 0.3e1_dp/ &
     585  1669325983 :                       0.2e1_dp*t131*rsrho + t21*t19*rsrho*t135
     586  1669325983 :                t139 = 0.1e1_dp/t26
     587  1669325983 :                t140 = t138*t139
     588  1669325983 :                e_c_u_0rho = -0.2e1_dp*t119*rsrho*t27 + t125*t140
     589  1669325983 :                t142 = t71**2
     590  1669325983 :                k_frho = t1/t142*t69/0.3e1_dp
     591  1669325983 :                t146 = 0.1e1_dp/t73
     592  1669325983 :                k_srho = t146*k_frho*t5
     593  1669325983 :                t148 = k_s**2
     594  1669325983 :                t149 = 0.1e1_dp/t148
     595  1669325983 :                t150 = my_norm_drho*t149
     596  1669325983 :                t151 = t6*k_srho
     597  1669325983 :                t153 = t75*t115
     598  1669325983 :                trho = -t150*t151/0.2e1_dp - t153/0.2e1_dp
     599  1669325983 :                t155 = gamma_var**2
     600  1669325983 :                t157 = beta/t155
     601  1669325983 :                t158 = t81**2
     602  1669325983 :                t159 = 0.1e1_dp/t158
     603  1669325983 :                Arho = t157*t159*e_c_u_0rho*t80
     604  1669325983 :                t162 = t78*t
     605  1669325983 :                t163 = t85*t91
     606  1669325983 :                t164 = t163*trho
     607  1669325983 :                t167 = Arho*t83
     608  1669325983 :                t168 = A*t
     609  1669325983 :                t170 = 0.2e1_dp*t168*trho
     610  1669325983 :                t171 = t167 + t170
     611  1669325983 :                t175 = t78*t83
     612  1669325983 :                t176 = t90**2
     613  1669325983 :                t177 = 0.1e1_dp/t176
     614  1669325983 :                t178 = t85*t177
     615  1669325983 :                t179 = A*t88
     616  1669325983 :                t182 = t83*t
     617  1669325983 :                t183 = t87*t182
     618  1669325983 :                t186 = t167 + t170 + 0.2e1_dp*t179*Arho + 0.4e1_dp*t183*trho
     619  1669325983 :                t187 = t178*t186
     620  1669325983 :                t189 = 0.2e1_dp*t162*t164 + t78*t83*t171*t91 - t175*t187
     621  1669325983 :                t190 = gamma_var*t189
     622  1669325983 :                t191 = 0.1e1_dp/t94
     623  1669325983 :                epsilon_cGGArho = e_c_u_0rho + t190*t191
     624  1669325983 :                kfrho = k_frho
     625  1669325983 :                ex_unifrho = -0.3e1_dp/0.4e1_dp*t5*kfrho
     626  1669325983 :                t195 = kf**2
     627  1669325983 :                t196 = 0.1e1_dp/t195
     628  1669325983 :                t197 = my_norm_drho*t196
     629  1669325983 :                t198 = t6*kfrho
     630  1669325983 :                t200 = t99*t115
     631  1669325983 :                srho = -t197*t198/0.2e1_dp - t200/0.2e1_dp
     632  1669325983 :                t202 = t105**2
     633  1669325983 :                t204 = 0.1e1_dp/t202*mu
     634  1669325983 :                Fxrho = 0.2e1_dp*t204*s*srho
     635  1669325983 :                t208 = my_rho*ex_unifrho
     636             : 
     637  1669325983 :                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  1575802362 :                               scale_ec*(epsilon_cGGA + my_rho*epsilon_cGGArho)
     641             :                END IF
     642             : 
     643  1669325983 :                tnorm_drho = t74*t6/0.2e1_dp
     644  1669325983 :                t214 = t163*tnorm_drho
     645  1669325983 :                t217 = t78*t182
     646  1669325983 :                t218 = A*tnorm_drho
     647  1669325983 :                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  1669325983 :                       t175*t178*t226
     650  1669325983 :                t230 = gamma_var*t229
     651  1669325983 :                Hnorm_drho = t230*t191
     652  1669325983 :                snorm_drho = t98*t6/0.2e1_dp
     653  1669325983 :                Fxnorm_drho = 0.2e1_dp*t204*s*snorm_drho
     654             : 
     655  1669325983 :                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  1575802362 :                                 Hnorm_drho
     659             :                END IF
     660             : 
     661  1669325983 :                t238 = 0.1e1_dp/t69
     662  1669325983 :                t239 = 0.1e1_dp/t111/t7*t238
     663  1669325983 :                t240 = t114**2
     664  1669325983 :                t241 = 0.1e1_dp/t240
     665  1669325983 :                t246 = 0.1e1_dp/t114/my_rho
     666             :                rsrhorho = -t4*t239*t241/0.18e2_dp + t4*t113* &
     667  1669325983 :                           t246/0.6e1_dp
     668  1669325983 :                t252 = 0.2e1_dp*t119*rsrhorho*t27
     669  1669325983 :                t253 = alpha_1_1*rsrho
     670  1669325983 :                t255 = t124*t138*t139
     671  1669325983 :                t259 = 0.1e1_dp/t123/t22
     672  1669325983 :                t260 = t11*t259
     673  1669325983 :                t261 = t138**2
     674  1669325983 :                t262 = t261*t139
     675  1669325983 :                t265 = 0.1e1_dp/t17
     676  1669325983 :                t266 = beta_1_1*t265
     677  1669325983 :                t267 = rsrho**2
     678  1669325983 :                t271 = t127*rsrhorho/0.2e1_dp
     679  1669325983 :                t272 = beta_2_1*rsrhorho
     680  1669325983 :                t273 = beta_3_1*t126
     681  1669325983 :                t277 = 0.3e1_dp/0.2e1_dp*t131*rsrhorho
     682  1669325983 :                t278 = t19**2
     683  1669325983 :                t280 = rs**2
     684  1669325983 :                t281 = 0.1e1_dp/t280
     685  1669325983 :                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  1669325983 :                       *t267*t281
     689  1669325983 :                t291 = t290*t139
     690  1669325983 :                t293 = t123**2
     691  1669325983 :                t294 = 0.1e1_dp/t293
     692  1669325983 :                t295 = t11*t294
     693  1669325983 :                t296 = t26**2
     694  1669325983 :                t297 = 0.1e1_dp/t296
     695  1669325983 :                t299 = t261*t297*t13
     696             :                e_c_u_0rhorho = -t252 + 0.2e1_dp*t253*t255 - 0.2e1_dp*t260* &
     697  1669325983 :                                t262 + t125*t291 + t295*t299/0.2e1_dp
     698  1669325983 :                e_c_u_01rho = e_c_u_0rho
     699  1669325983 :                t305 = t69**2
     700  1669325983 :                k_frhorho = -0.2e1_dp/0.9e1_dp*t1/t142/t70*t305
     701  1669325983 :                t309 = 0.1e1_dp/t73/t72
     702  1669325983 :                t310 = k_frho**2
     703  1669325983 :                t315 = t146*k_frhorho*t5
     704  1669325983 :                k_srhorho = -t309*t310*t238/0.2e1_dp + t315
     705  1669325983 :                k_s1rho = k_srho
     706  1669325983 :                t317 = 0.1e1_dp/t148/k_s
     707  1669325983 :                t318 = my_norm_drho*t317
     708  1669325983 :                t321 = t115*k_srho
     709  1669325983 :                t323 = t150*t321/0.2e1_dp
     710  1669325983 :                t324 = t6*k_srhorho
     711  1669325983 :                t327 = t115*k_s1rho
     712  1669325983 :                t329 = t150*t327/0.2e1_dp
     713  1669325983 :                t330 = t75*t246
     714             :                trhorho = t318*t151*k_s1rho + t323 - t150*t324/0.2e1_dp + &
     715  1669325983 :                          t329 + t330
     716  1669325983 :                t331 = t6*k_s1rho
     717  1669325983 :                t1rho = -t150*t331/0.2e1_dp - t153/0.2e1_dp
     718  1669325983 :                t336 = beta/t155/gamma_var
     719  1669325983 :                t338 = 0.1e1_dp/t158/t81
     720  1669325983 :                t339 = t336*t338
     721  1669325983 :                t340 = t80**2
     722  1669325983 :                t341 = e_c_u_0rho*t340
     723  1669325983 :                t348 = t336*t159
     724  1669325983 :                t349 = e_c_u_0rho*e_c_u_01rho
     725             :                Arhorho = 0.2e1_dp*t339*t341*e_c_u_01rho + t157*t159* &
     726  1669325983 :                          e_c_u_0rhorho*t80 - t348*t349*t80
     727  1669325983 :                A1rho = t157*t159*e_c_u_01rho*t80
     728  1669325983 :                t354 = t78*t1rho
     729  1669325983 :                t357 = A1rho*t83
     730  1669325983 :                t359 = 0.2e1_dp*t168*t1rho
     731  1669325983 :                t360 = t357 + t359
     732  1669325983 :                t361 = t360*t91
     733  1669325983 :                t362 = t361*trho
     734  1669325983 :                t369 = t357 + t359 + 0.2e1_dp*t179*A1rho + 0.4e1_dp*t183*t1rho
     735  1669325983 :                t370 = trho*t369
     736  1669325983 :                t371 = t178*t370
     737  1669325983 :                t374 = t163*trhorho
     738  1669325983 :                t377 = t171*t91
     739  1669325983 :                t378 = t377*t1rho
     740  1669325983 :                t381 = Arhorho*t83
     741  1669325983 :                t382 = Arho*t
     742  1669325983 :                t384 = 0.2e1_dp*t382*t1rho
     743  1669325983 :                t385 = A1rho*t
     744  1669325983 :                t387 = 0.2e1_dp*t385*trho
     745  1669325983 :                t388 = A*t1rho
     746  1669325983 :                t390 = 0.2e1_dp*t388*trho
     747  1669325983 :                t392 = 0.2e1_dp*t168*trhorho
     748  1669325983 :                t393 = t381 + t384 + t387 + t390 + t392
     749  1669325983 :                t397 = t171*t177
     750  1669325983 :                t400 = t186*t1rho
     751  1669325983 :                t401 = t178*t400
     752  1669325983 :                t404 = t360*t177
     753  1669325983 :                t408 = 0.1e1_dp/t176/t90
     754  1669325983 :                t409 = t85*t408
     755  1669325983 :                t410 = t186*t369
     756  1669325983 :                t414 = A1rho*t88
     757  1669325983 :                t417 = A*t182
     758  1669325983 :                t418 = Arho*t1rho
     759  1669325983 :                t423 = trho*A1rho
     760  1669325983 :                t426 = t87*t83
     761  1669325983 :                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  1669325983 :                       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  1669325983 :                       *t432
     770  1669325983 :                t436 = gamma_var*t435
     771  1669325983 :                t438 = t94**2
     772  1669325983 :                t439 = 0.1e1_dp/t438
     773  1669325983 :                t440 = t163*t1rho
     774             :                t448 = 0.2e1_dp*t162*t440 + t78*t83*t360*t91 - t175* &
     775  1669325983 :                       t178*t369
     776  1669325983 :                t449 = t439*t448
     777  1669325983 :                t451 = gamma_var*t448
     778  1669325983 :                epsilon_cGGArhorho = e_c_u_0rhorho + t436*t191 - t190*t449
     779  1669325983 :                kfrhorho = k_frhorho
     780  1669325983 :                ex_unifrhorho = -0.3e1_dp/0.4e1_dp*t5*kfrhorho
     781  1669325983 :                ex_unif1rho = ex_unifrho
     782  1669325983 :                t456 = 0.1e1_dp/t195/kf
     783  1669325983 :                t457 = my_norm_drho*t456
     784  1669325983 :                t458 = kfrho**2
     785  1669325983 :                t459 = t6*t458
     786  1669325983 :                t461 = t115*kfrho
     787  1669325983 :                t462 = t197*t461
     788  1669325983 :                t463 = t6*kfrhorho
     789  1669325983 :                t465 = t197*t463/0.2e1_dp
     790  1669325983 :                t466 = t99*t246
     791  1669325983 :                srhorho = t457*t459 + t462 - t465 + t466
     792  1669325983 :                s1rho = srho
     793  1669325983 :                t469 = mu**2
     794  1669325983 :                t470 = 0.1e1_dp/t202/t105*t469
     795  1669325983 :                t471 = t470*t101
     796  1669325983 :                t472 = srho*t103
     797  1669325983 :                t476 = s1rho*srho
     798             :                Fxrhorho = -0.8e1_dp*t471*t472*s1rho + 0.2e1_dp*t204* &
     799  1669325983 :                           t476 + 0.2e1_dp*t204*s*srhorho
     800  1669325983 :                Fx1rho = 0.2e1_dp*t204*s*s1rho
     801  1669325983 :                t487 = my_rho*ex_unifrhorho
     802  1669325983 :                t491 = my_rho*ex_unif1rho
     803             : 
     804  1669325983 :                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   185300121 :                                                                                 + epsilon_cGGArho + my_rho*epsilon_cGGArhorho)
     810             :                END IF
     811             : 
     812  1669325983 :                t496 = t149*t6
     813  1669325983 :                t498 = t74*t115
     814  1669325983 :                tnorm_drhorho = -t496*k_srho/0.2e1_dp - t498/0.2e1_dp
     815  1669325983 :                t500 = t78*trho
     816  1669325983 :                t503 = t377*tnorm_drho
     817  1669325983 :                t506 = tnorm_drho*t186
     818  1669325983 :                t507 = t178*t506
     819  1669325983 :                t510 = t163*tnorm_drhorho
     820  1669325983 :                t513 = t91*trho
     821  1669325983 :                t517 = Arho*tnorm_drho
     822  1669325983 :                t521 = A*tnorm_drhorho
     823  1669325983 :                t525 = t177*t186
     824  1669325983 :                t529 = t226*trho
     825  1669325983 :                t530 = t178*t529
     826  1669325983 :                t535 = t226*t186
     827  1669325983 :                t541 = A*trho
     828  1669325983 :                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  1669325983 :                       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  1669325983 :                       t397*t226 + 0.2e1_dp*t175*t409*t535 - t175*t178*t553
     837  1669325983 :                t557 = gamma_var*t556
     838  1669325983 :                t559 = t439*t189
     839  1669325983 :                Hnorm_drhorho = t557*t191 - t230*t559
     840  1669325983 :                t562 = t196*t6
     841  1669325983 :                snorm_drhorho = -t562*kfrho/0.2e1_dp - t98*t115/0.2e1_dp
     842  1669325983 :                t566 = snorm_drho*t103
     843             :                Fxnorm_drhorho = -0.8e1_dp*t471*t566*srho + 0.2e1_dp*t204 &
     844  1669325983 :                                 *srho*snorm_drho + 0.2e1_dp*t204*s*snorm_drhorho
     845             : 
     846  1669325983 :                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   185300121 :                                                                                              *Hnorm_drhorho)
     851             :                END IF
     852             : 
     853  1669325983 :                t581 = tnorm_drho**2
     854  1669325983 :                t586 = A*t581
     855  1669325983 :                t590 = tnorm_drho*t226
     856  1669325983 :                t591 = t178*t590
     857  1669325983 :                t594 = t177*t226
     858  1669325983 :                t598 = t226**2
     859  1669325983 :                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  1669325983 :                                  t218*t594 + 0.2e1_dp*t175*t409*t598 - t175*t178*t605)
     863  1669325983 :                t611 = t229**2
     864  1669325983 :                t612 = gamma_var*t611
     865  1669325983 :                Hnorm_drhonorm_drho = t609*t191 - t612*t439
     866  1669325983 :                t614 = snorm_drho**2
     867             :                Fxnorm_drhonorm_drho = -0.8e1_dp*t470*t101*t614*t103 + &
     868  1669325983 :                                       0.2e1_dp*t204*t614
     869             : 
     870  1669325983 :                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   185300121 :                                       scale_ec*my_rho*Hnorm_drhonorm_drho
     874             :                END IF
     875             : 
     876  1669325983 :                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       84378 :    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       90125 :    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       18025 :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, &
    1337       18025 :          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       18025 :          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       18025 :       CALL timeset(routineN, handle)
    1342       18025 :       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       18025 :                           local_bounds=bo)
    1349       18025 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
    1350             : 
    1351       18025 :       dummy => rhoa
    1352       18025 :       e_0 => dummy
    1353       18025 :       e_ra => dummy
    1354       18025 :       e_rb => dummy
    1355       18025 :       e_ndra_ra => dummy
    1356       18025 :       e_ndrb_rb => dummy
    1357       18025 :       e_ndr_ndr => dummy
    1358       18025 :       e_ndra_ndra => dummy
    1359       18025 :       e_ndrb_ndrb => dummy
    1360       18025 :       e_ndr => dummy
    1361       18025 :       e_ndra => dummy
    1362       18025 :       e_ndrb => dummy
    1363       18025 :       e_ra_ra => dummy
    1364       18025 :       e_ra_rb => dummy
    1365       18025 :       e_rb_rb => dummy
    1366       18025 :       e_ndr_ra => dummy
    1367       18025 :       e_ndr_rb => dummy
    1368             : 
    1369       18025 :       IF (grad_deriv >= 0) THEN
    1370             :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
    1371       18025 :                                          allocate_deriv=.TRUE.)
    1372       18025 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
    1373             :       END IF
    1374       18025 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1375             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
    1376       16151 :                                          allocate_deriv=.TRUE.)
    1377       16151 :          CALL xc_derivative_get(deriv, deriv_data=e_ra)
    1378             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
    1379       16151 :                                          allocate_deriv=.TRUE.)
    1380       16151 :          CALL xc_derivative_get(deriv, deriv_data=e_rb)
    1381             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
    1382       16151 :                                          allocate_deriv=.TRUE.)
    1383       16151 :          CALL xc_derivative_get(deriv, deriv_data=e_ndr)
    1384             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
    1385       16151 :                                          allocate_deriv=.TRUE.)
    1386       16151 :          CALL xc_derivative_get(deriv, deriv_data=e_ndra)
    1387             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
    1388       16151 :                                          allocate_deriv=.TRUE.)
    1389       16151 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
    1390             :       END IF
    1391       18025 :       IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
    1392             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
    1393        1214 :                                          allocate_deriv=.TRUE.)
    1394        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_ra_ra)
    1395             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
    1396        1214 :                                          allocate_deriv=.TRUE.)
    1397        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_ra_rb)
    1398             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
    1399        1214 :                                          allocate_deriv=.TRUE.)
    1400        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_rb_rb)
    1401             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
    1402        1214 :                                          allocate_deriv=.TRUE.)
    1403        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_ndr_ra)
    1404             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
    1405        1214 :                                          allocate_deriv=.TRUE.)
    1406        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_ndr_rb)
    1407             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
    1408        1214 :                                          allocate_deriv=.TRUE.)
    1409        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_ndra_ra)
    1410             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
    1411        1214 :                                          allocate_deriv=.TRUE.)
    1412        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrb_rb)
    1413             :          deriv => xc_dset_get_derivative(deriv_set, &
    1414        1214 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
    1415        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_ndr_ndr)
    1416             :          deriv => xc_dset_get_derivative(deriv_set, &
    1417        1214 :                                          [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
    1418        1214 :          CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndra)
    1419             :          deriv => xc_dset_get_derivative(deriv_set, &
    1420        1214 :                                          [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
    1421        1214 :          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       18025 :       CALL section_vals_val_get(pbe_params, "scale_c", r_val=scale_ec)
    1428       18025 :       CALL section_vals_val_get(pbe_params, "scale_x", r_val=scale_ex)
    1429       18025 :       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       18025 : !$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       18025 :       CALL timestop(handle)
    1451       18025 :    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       18025 :    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       18025 :       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       18025 :          CPABORT("Unsupported parametrization")
    1573             :       END SELECT
    1574             : 
    1575       18025 :       gamma_var = (0.1e1_dp - LOG(0.2e1_dp))/pi**2
    1576       18025 :       p_1 = 0.10e1_dp
    1577       18025 :       A_1 = 0.31091e-1_dp
    1578       18025 :       alpha_1_1 = 0.21370e0_dp
    1579       18025 :       beta_1_1 = 0.75957e1_dp
    1580       18025 :       beta_2_1 = 0.35876e1_dp
    1581       18025 :       beta_3_1 = 0.16382e1_dp
    1582       18025 :       beta_4_1 = 0.49294e0_dp
    1583       18025 :       p_2 = 0.10e1_dp
    1584       18025 :       A_2 = 0.15545e-1_dp
    1585       18025 :       alpha_1_2 = 0.20548e0_dp
    1586       18025 :       beta_1_2 = 0.141189e2_dp
    1587       18025 :       beta_2_2 = 0.61977e1_dp
    1588       18025 :       beta_3_2 = 0.33662e1_dp
    1589       18025 :       beta_4_2 = 0.62517e0_dp
    1590       18025 :       p_3 = 0.10e1_dp
    1591       18025 :       A_3 = 0.16887e-1_dp
    1592       18025 :       alpha_1_3 = 0.11125e0_dp
    1593       18025 :       beta_1_3 = 0.10357e2_dp
    1594       18025 :       beta_2_3 = 0.36231e1_dp
    1595       18025 :       beta_3_3 = 0.88026e0_dp
    1596       18025 :       beta_4_3 = 0.49671e0_dp
    1597       18025 :       t3 = 3**(0.1e1_dp/0.3e1_dp)
    1598       18025 :       t4 = 4**(0.1e1_dp/0.3e1_dp)
    1599       18025 :       t5 = t4**2
    1600       18025 :       t6 = t3*t5
    1601       18025 :       t7 = 0.1e1_dp/pi
    1602       18025 :       t90 = pi**2
    1603             : 
    1604       18025 :       SELECT CASE (grad_deriv)
    1605             :       CASE default
    1606       18025 : !$OMP DO
    1607             :          DO ii = 1, npoints
    1608   580893420 :             my_rhoa = MAX(rhoa(ii), 0.0_dp)
    1609   580893420 :             my_rhob = MAX(rhob(ii), 0.0_dp)
    1610   580893420 :             my_rho = my_rhoa + my_rhob
    1611   580893420 :             IF (my_rho > epsilon_rho) THEN
    1612   534627218 :                my_rhoa = MAX(EPSILON(0.0_dp)*1.e4_dp, my_rhoa)
    1613   534627218 :                my_rhob = MAX(EPSILON(0.0_dp)*1.e4_dp, my_rhob)
    1614   534627218 :                my_rho = my_rhoa + my_rhob
    1615   534627218 :                my_norm_drho = norm_drho(ii)
    1616   534627218 :                my_norm_drhoa = norm_drhoa(ii)
    1617   534627218 :                my_norm_drhob = norm_drhob(ii)
    1618             : 
    1619   534627218 :                t1 = my_rhoa - my_rhob
    1620   534627218 :                t2 = 0.1e1_dp/my_rho
    1621   534627218 :                chi = t1*t2
    1622   534627218 :                t8 = t7*t2
    1623   534627218 :                t9 = t8**(0.1e1_dp/0.3e1_dp)
    1624   534627218 :                rs = t6*t9/0.4e1_dp
    1625   534627218 :                t12 = 0.1e1_dp + alpha_1_1*rs
    1626   534627218 :                t14 = 0.1e1_dp/A_1
    1627   534627218 :                t15 = SQRT(rs)
    1628   534627218 :                t18 = t15*rs
    1629   534627218 :                t20 = p_1 + 0.1e1_dp
    1630   534627218 :                t21 = rs**t20
    1631   534627218 :                t22 = beta_4_1*t21
    1632   534627218 :                t23 = beta_1_1*t15 + beta_2_1*rs + beta_3_1*t18 + t22
    1633   534627218 :                t27 = 0.1e1_dp + t14/t23/0.2e1_dp
    1634   534627218 :                t28 = LOG(t27)
    1635   534627218 :                e_c_u_0 = -0.2e1_dp*A_1*t12*t28
    1636   534627218 :                t32 = 0.1e1_dp + alpha_1_2*rs
    1637   534627218 :                t34 = 0.1e1_dp/A_2
    1638   534627218 :                t38 = p_2 + 0.1e1_dp
    1639   534627218 :                t39 = rs**t38
    1640   534627218 :                t40 = beta_4_2*t39
    1641   534627218 :                t41 = beta_1_2*t15 + beta_2_2*rs + beta_3_2*t18 + t40
    1642   534627218 :                t45 = 0.1e1_dp + t34/t41/0.2e1_dp
    1643   534627218 :                t46 = LOG(t45)
    1644   534627218 :                t50 = 0.1e1_dp + alpha_1_3*rs
    1645   534627218 :                t52 = 0.1e1_dp/A_3
    1646   534627218 :                t56 = p_3 + 0.1e1_dp
    1647   534627218 :                t57 = rs**t56
    1648   534627218 :                t58 = beta_4_3*t57
    1649   534627218 :                t59 = beta_1_3*t15 + beta_2_3*rs + beta_3_3*t18 + t58
    1650   534627218 :                t63 = 0.1e1_dp + t52/t59/0.2e1_dp
    1651   534627218 :                t64 = LOG(t63)
    1652   534627218 :                alpha_c = 0.2e1_dp*A_3*t50*t64
    1653   534627218 :                t66 = 2**(0.1e1_dp/0.3e1_dp)
    1654   534627218 :                t69 = 1/(2*t66 - 2)
    1655   534627218 :                t70 = 0.1e1_dp + chi
    1656   534627218 :                t71 = t70**(0.1e1_dp/0.3e1_dp)
    1657   534627218 :                t72 = t71*t70
    1658   534627218 :                t73 = 0.1e1_dp - chi
    1659   534627218 :                t74 = t73**(0.1e1_dp/0.3e1_dp)
    1660   534627218 :                t75 = t74*t73
    1661   534627218 :                f = (t72 + t75 - 0.2e1_dp)*t69
    1662   534627218 :                t77 = alpha_c*f
    1663   534627218 :                t78 = 0.9e1_dp/0.8e1_dp/t69
    1664   534627218 :                t79 = chi**2
    1665   534627218 :                t80 = t79**2
    1666   534627218 :                t82 = t78*(0.1e1_dp - t80)
    1667   534627218 :                t84 = -0.2e1_dp*A_2*t32*t46 - e_c_u_0
    1668   534627218 :                t85 = t84*f
    1669   534627218 :                epsilon_c_unif = e_c_u_0 + t77*t82 + t85*t80
    1670   534627218 :                t87 = t71**2
    1671   534627218 :                t88 = t74**2
    1672   534627218 :                phi = t87/0.2e1_dp + t88/0.2e1_dp
    1673   534627218 :                t91 = t90*my_rho
    1674   534627218 :                t92 = t91**(0.1e1_dp/0.3e1_dp)
    1675   534627218 :                t93 = t3*t92*t7
    1676   534627218 :                t94 = SQRT(t93)
    1677   534627218 :                k_s = 0.2e1_dp*t94
    1678   534627218 :                t95 = 0.1e1_dp/phi
    1679   534627218 :                t96 = my_norm_drho*t95
    1680   534627218 :                t97 = 0.1e1_dp/k_s
    1681   534627218 :                t98 = t97*t2
    1682   534627218 :                t = t96*t98/0.2e1_dp
    1683   534627218 :                t100 = 0.1e1_dp/gamma_var
    1684   534627218 :                t101 = beta*t100
    1685   534627218 :                t102 = epsilon_c_unif*t100
    1686   534627218 :                t103 = phi**2
    1687   534627218 :                t104 = t103*phi
    1688   534627218 :                t105 = 0.1e1_dp/t104
    1689   534627218 :                t107 = EXP(-t102*t105)
    1690   534627218 :                t108 = t107 - 0.1e1_dp
    1691   534627218 :                A = t101/t108
    1692   534627218 :                t110 = gamma_var*t104
    1693   534627218 :                t111 = t**2
    1694   534627218 :                t112 = A*t111
    1695   534627218 :                t113 = 0.1e1_dp + t112
    1696   534627218 :                t115 = A**2
    1697   534627218 :                t116 = t111**2
    1698   534627218 :                t118 = 0.1e1_dp + t112 + t115*t116
    1699   534627218 :                t119 = 0.1e1_dp/t118
    1700   534627218 :                t122 = 0.1e1_dp + t101*t111*t113*t119
    1701   534627218 :                t123 = LOG(t122)
    1702   534627218 :                epsilon_cGGA = epsilon_c_unif + t110*t123
    1703   534627218 :                t124 = t3*t66
    1704   534627218 :                t125 = t90*my_rhoa
    1705   534627218 :                t126 = t125**(0.1e1_dp/0.3e1_dp)
    1706   534627218 :                kf_a = t124*t126
    1707   534627218 :                ex_unif_a = -0.3e1_dp/0.4e1_dp*t7*kf_a
    1708   534627218 :                t129 = 0.1e1_dp/kf_a
    1709   534627218 :                t130 = my_norm_drhoa*t129
    1710   534627218 :                t131 = 0.1e1_dp/my_rhoa
    1711   534627218 :                s_a = t130*t131/0.2e1_dp
    1712   534627218 :                t133 = s_a**2
    1713   534627218 :                t135 = 0.1e1_dp/kappa
    1714   534627218 :                t137 = 0.1e1_dp + mu*t133*t135
    1715   534627218 :                Fx_a = 0.1e1_dp + kappa - kappa/t137
    1716   534627218 :                t140 = my_rhoa*ex_unif_a
    1717   534627218 :                t142 = t90*my_rhob
    1718   534627218 :                t143 = t142**(0.1e1_dp/0.3e1_dp)
    1719   534627218 :                kf_b = t124*t143
    1720   534627218 :                ex_unif_b = -0.3e1_dp/0.4e1_dp*t7*kf_b
    1721   534627218 :                t146 = 0.1e1_dp/kf_b
    1722   534627218 :                t147 = my_norm_drhob*t146
    1723   534627218 :                t148 = 0.1e1_dp/my_rhob
    1724   534627218 :                s_b = t147*t148/0.2e1_dp
    1725   534627218 :                t150 = s_b**2
    1726   534627218 :                t153 = 0.1e1_dp + mu*t150*t135
    1727   534627218 :                Fx_b = 0.1e1_dp + kappa - kappa/t153
    1728   534627218 :                t156 = my_rhob*ex_unif_b
    1729             : 
    1730   534627218 :                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   534627218 :                             /0.2e1_dp + scale_ec*my_rho*epsilon_cGGA
    1734             :                END IF
    1735             : 
    1736   534627218 :                t162 = my_rho**2
    1737   534627218 :                t163 = 0.1e1_dp/t162
    1738   534627218 :                t164 = t1*t163
    1739   534627218 :                chirhoa = t2 - t164
    1740   534627218 :                t165 = t9**2
    1741   534627218 :                t167 = 0.1e1_dp/t165*t7
    1742   534627218 :                rsrhoa = -t6*t167*t163/0.12e2_dp
    1743   534627218 :                t171 = A_1*alpha_1_1
    1744   534627218 :                t175 = t23**2
    1745   534627218 :                t176 = 0.1e1_dp/t175
    1746   534627218 :                t177 = t12*t176
    1747   534627218 :                t178 = 0.1e1_dp/t15
    1748   534627218 :                t179 = beta_1_1*t178
    1749   534627218 :                t183 = beta_3_1*t15
    1750   534627218 :                t187 = 0.1e1_dp/rs
    1751             :                t190 = t179*rsrhoa/0.2e1_dp + beta_2_1*rsrhoa + 0.3e1_dp/ &
    1752   534627218 :                       0.2e1_dp*t183*rsrhoa + t22*t20*rsrhoa*t187
    1753   534627218 :                t191 = 0.1e1_dp/t27
    1754   534627218 :                t192 = t190*t191
    1755   534627218 :                e_c_u_0rhoa = -0.2e1_dp*t171*rsrhoa*t28 + t177*t192
    1756   534627218 :                t194 = A_2*alpha_1_2
    1757   534627218 :                t198 = t41**2
    1758   534627218 :                t199 = 0.1e1_dp/t198
    1759   534627218 :                t200 = t32*t199
    1760   534627218 :                t201 = beta_1_2*t178
    1761   534627218 :                t205 = beta_3_2*t15
    1762             :                t211 = t201*rsrhoa/0.2e1_dp + beta_2_2*rsrhoa + 0.3e1_dp/ &
    1763   534627218 :                       0.2e1_dp*t205*rsrhoa + t40*t38*rsrhoa*t187
    1764   534627218 :                t212 = 0.1e1_dp/t45
    1765   534627218 :                t213 = t211*t212
    1766   534627218 :                e_c_u_1rhoa = -0.2e1_dp*t194*rsrhoa*t46 + t200*t213
    1767   534627218 :                t215 = A_3*alpha_1_3
    1768   534627218 :                t219 = t59**2
    1769   534627218 :                t220 = 0.1e1_dp/t219
    1770   534627218 :                t221 = t50*t220
    1771   534627218 :                t222 = beta_1_3*t178
    1772   534627218 :                t226 = beta_3_3*t15
    1773             :                t232 = t222*rsrhoa/0.2e1_dp + beta_2_3*rsrhoa + 0.3e1_dp/ &
    1774   534627218 :                       0.2e1_dp*t226*rsrhoa + t58*t56*rsrhoa*t187
    1775   534627218 :                t233 = 0.1e1_dp/t63
    1776   534627218 :                t234 = t232*t233
    1777   534627218 :                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   534627218 :                         *t74*chirhoa)*t69
    1780   534627218 :                t240 = alpha_crhoa*f
    1781   534627218 :                t242 = alpha_c*frhoa
    1782   534627218 :                t244 = t79*chi
    1783   534627218 :                t245 = t78*t244
    1784   534627218 :                t246 = t245*chirhoa
    1785   534627218 :                t248 = 0.4e1_dp*t77*t246
    1786   534627218 :                t249 = e_c_u_1rhoa - e_c_u_0rhoa
    1787   534627218 :                t250 = t249*f
    1788   534627218 :                t252 = t84*frhoa
    1789   534627218 :                t254 = t244*chirhoa
    1790   534627218 :                t256 = 0.4e1_dp*t85*t254
    1791             :                epsilon_c_unifrhoa = e_c_u_0rhoa + t240*t82 + t242*t82 - t248 &
    1792   534627218 :                                     + t250*t80 + t252*t80 + t256
    1793   534627218 :                t257 = 0.1e1_dp/t71
    1794   534627218 :                t259 = 0.1e1_dp/t74
    1795   534627218 :                phirhoa = t257*chirhoa/0.3e1_dp - t259*chirhoa/0.3e1_dp
    1796   534627218 :                t262 = t92**2
    1797   534627218 :                k_frhoa = t3/t262*t90/0.3e1_dp
    1798   534627218 :                t266 = 0.1e1_dp/t94
    1799   534627218 :                k_srhoa = t266*k_frhoa*t7
    1800   534627218 :                t268 = 0.1e1_dp/t103
    1801   534627218 :                t269 = my_norm_drho*t268
    1802   534627218 :                t272 = k_s**2
    1803   534627218 :                t273 = 0.1e1_dp/t272
    1804   534627218 :                t274 = t273*t2
    1805   534627218 :                t277 = t97*t163
    1806   534627218 :                t278 = t96*t277
    1807             :                trhoa = -t269*t98*phirhoa/0.2e1_dp - t96*t274*k_srhoa/ &
    1808   534627218 :                        0.2e1_dp - t278/0.2e1_dp
    1809   534627218 :                t280 = t108**2
    1810   534627218 :                t281 = 0.1e1_dp/t280
    1811   534627218 :                t282 = epsilon_c_unifrhoa*t100
    1812   534627218 :                t284 = t103**2
    1813   534627218 :                t285 = 0.1e1_dp/t284
    1814   534627218 :                t286 = t285*phirhoa
    1815   534627218 :                t289 = -t282*t105 + 0.3e1_dp*t102*t286
    1816   534627218 :                Arhoa = -t101*t281*t289*t107
    1817   534627218 :                t293 = gamma_var*t103
    1818   534627218 :                t294 = t123*phirhoa
    1819   534627218 :                t297 = t101*t
    1820   534627218 :                t298 = t113*t119
    1821   534627218 :                t299 = t298*trhoa
    1822   534627218 :                t302 = Arhoa*t111
    1823   534627218 :                t303 = A*t
    1824   534627218 :                t305 = 0.2e1_dp*t303*trhoa
    1825   534627218 :                t306 = t302 + t305
    1826   534627218 :                t310 = t101*t111
    1827   534627218 :                t311 = t118**2
    1828   534627218 :                t312 = 0.1e1_dp/t311
    1829   534627218 :                t313 = t113*t312
    1830   534627218 :                t314 = A*t116
    1831   534627218 :                t317 = t111*t
    1832   534627218 :                t318 = t115*t317
    1833   534627218 :                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   534627218 :                       t313*t321
    1836   534627218 :                t325 = 0.1e1_dp/t122
    1837   534627218 :                t326 = t324*t325
    1838             :                epsilon_cGGArhoa = epsilon_c_unifrhoa + 0.3e1_dp*t293*t294 + &
    1839   534627218 :                                   t110*t326
    1840   534627218 :                t329 = t126**2
    1841   534627218 :                kf_arhoa = t124/t329*t90/0.3e1_dp
    1842   534627218 :                ex_unif_arhoa = -0.3e1_dp/0.4e1_dp*t7*kf_arhoa
    1843   534627218 :                t335 = kf_a**2
    1844   534627218 :                t336 = 0.1e1_dp/t335
    1845   534627218 :                t337 = my_norm_drhoa*t336
    1846   534627218 :                t340 = my_rhoa**2
    1847   534627218 :                t341 = 0.1e1_dp/t340
    1848   534627218 :                s_arhoa = -t337*t131*kf_arhoa/0.2e1_dp - t130*t341/0.2e1_dp
    1849   534627218 :                t344 = t137**2
    1850   534627218 :                t346 = 0.1e1_dp/t344*mu
    1851   534627218 :                Fx_arhoa = 0.2e1_dp*t346*s_a*s_arhoa
    1852   534627218 :                t350 = my_rhoa*ex_unif_arhoa
    1853             : 
    1854   534627218 :                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   454725488 :                              epsilon_cGGA + my_rho*epsilon_cGGArhoa)
    1859             :                END IF
    1860             : 
    1861   534627218 :                chirhob = -t2 - t164
    1862   534627218 :                rsrhob = rsrhoa
    1863             :                t368 = t179*rsrhob/0.2e1_dp + beta_2_1*rsrhob + 0.3e1_dp/ &
    1864   534627218 :                       0.2e1_dp*t183*rsrhob + t22*t20*rsrhob*t187
    1865   534627218 :                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   534627218 :                       0.2e1_dp*t205*rsrhob + t40*t38*rsrhob*t187
    1868   534627218 :                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   534627218 :                       0.2e1_dp*t226*rsrhob + t58*t56*rsrhob*t187
    1871   534627218 :                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   534627218 :                         *t74*chirhob)*t69
    1874   534627218 :                t403 = alpha_crhob*f
    1875   534627218 :                t405 = alpha_c*frhob
    1876   534627218 :                t407 = t245*chirhob
    1877   534627218 :                t409 = 0.4e1_dp*t77*t407
    1878   534627218 :                t410 = e_c_u_1rhob - e_c_u_0rhob
    1879   534627218 :                t411 = t410*f
    1880   534627218 :                t413 = t84*frhob
    1881   534627218 :                t415 = t244*chirhob
    1882   534627218 :                t417 = 0.4e1_dp*t85*t415
    1883             :                epsilon_c_unifrhob = e_c_u_0rhob + t403*t82 + t405*t82 - t409 &
    1884   534627218 :                                     + t411*t80 + t413*t80 + t417
    1885   534627218 :                phirhob = t257*chirhob/0.3e1_dp - t259*chirhob/0.3e1_dp
    1886   534627218 :                k_frhob = k_frhoa
    1887   534627218 :                k_srhob = t266*k_frhob*t7
    1888             :                trhob = -t269*t98*phirhob/0.2e1_dp - t96*t274*k_srhob/ &
    1889   534627218 :                        0.2e1_dp - t278/0.2e1_dp
    1890   534627218 :                t427 = epsilon_c_unifrhob*t100
    1891   534627218 :                t429 = t285*phirhob
    1892   534627218 :                t432 = -t427*t105 + 0.3e1_dp*t102*t429
    1893   534627218 :                Arhob = -t101*t281*t432*t107
    1894   534627218 :                t436 = t123*phirhob
    1895   534627218 :                t439 = t298*trhob
    1896   534627218 :                t442 = Arhob*t111
    1897   534627218 :                t444 = 0.2e1_dp*t303*trhob
    1898   534627218 :                t445 = t442 + t444
    1899   534627218 :                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   534627218 :                       t313*t453
    1902   534627218 :                t457 = t456*t325
    1903             :                epsilon_cGGArhob = epsilon_c_unifrhob + 0.3e1_dp*t293*t436 + &
    1904   534627218 :                                   t110*t457
    1905   534627218 :                t460 = t143**2
    1906   534627218 :                kf_brhob = t124/t460*t90/0.3e1_dp
    1907   534627218 :                ex_unif_brhob = -0.3e1_dp/0.4e1_dp*t7*kf_brhob
    1908   534627218 :                t466 = kf_b**2
    1909   534627218 :                t467 = 0.1e1_dp/t466
    1910   534627218 :                t468 = my_norm_drhob*t467
    1911   534627218 :                t471 = my_rhob**2
    1912   534627218 :                t472 = 0.1e1_dp/t471
    1913   534627218 :                s_brhob = -t468*t148*kf_brhob/0.2e1_dp - t147*t472/0.2e1_dp
    1914   534627218 :                t475 = t153**2
    1915   534627218 :                t477 = 0.1e1_dp/t475*mu
    1916   534627218 :                Fx_brhob = 0.2e1_dp*t477*s_b*s_brhob
    1917   534627218 :                t481 = my_rhob*ex_unif_brhob
    1918             : 
    1919   534627218 :                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   454725488 :                              epsilon_cGGA + my_rho*epsilon_cGGArhob)
    1924             :                END IF
    1925             : 
    1926   534627218 :                t488 = t95*t97
    1927   534627218 :                tnorm_drho = t488*t2/0.2e1_dp
    1928   534627218 :                t493 = t101*t317
    1929   534627218 :                t494 = A*tnorm_drho
    1930   534627218 :                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   534627218 :                       t494*t119 - t310*t313*t502
    1933   534627218 :                t506 = t505*t325
    1934   534627218 :                Hnorm_drho = t110*t506
    1935             : 
    1936   534627218 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1937             :                   e_ndr(ii) = e_ndr(ii) + &
    1938   454725488 :                               scale_ec*my_rho*Hnorm_drho
    1939             :                END IF
    1940             : 
    1941   534627218 :                s_anorm_drhoa = t129*t131/0.2e1_dp
    1942   534627218 :                Fx_anorm_drhoa = 0.2e1_dp*t346*s_a*s_anorm_drhoa
    1943             : 
    1944   534627218 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1945             :                   e_ndra(ii) = e_ndra(ii) + &
    1946   454725488 :                                scale_ex*t140*Fx_anorm_drhoa
    1947             :                END IF
    1948             : 
    1949   534627218 :                s_bnorm_drhob = t146*t148/0.2e1_dp
    1950   534627218 :                Fx_bnorm_drhob = 0.2e1_dp*t477*s_b*s_bnorm_drhob
    1951             : 
    1952   534627218 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
    1953             :                   e_ndrb(ii) = e_ndrb(ii) + &
    1954   454725488 :                                scale_ex*t156*Fx_bnorm_drhob
    1955             :                END IF
    1956             : 
    1957   534627218 :                IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
    1958    18956493 :                   t518 = 0.1e1_dp/t162/my_rho
    1959    18956493 :                   t519 = t1*t518
    1960    18956493 :                   chirhoarhoa = -0.2e1_dp*t163 + 0.2e1_dp*t519
    1961    18956493 :                   t523 = 0.1e1_dp/t90
    1962    18956493 :                   t525 = t162**2
    1963             :                   rsrhoarhoa = -t6/t165/t8*t523/t525/0.18e2_dp + &
    1964    18956493 :                                t6*t167*t518/0.6e1_dp
    1965    18956493 :                   t536 = alpha_1_1*rsrhoa
    1966    18956493 :                   t538 = t176*t190*t191
    1967    18956493 :                   t543 = t12/t175/t23
    1968    18956493 :                   t544 = t190**2
    1969    18956493 :                   t548 = 0.1e1_dp/t18
    1970    18956493 :                   t549 = beta_1_1*t548
    1971    18956493 :                   t550 = rsrhoa**2
    1972    18956493 :                   t556 = beta_3_1*t178
    1973    18956493 :                   t561 = t20**2
    1974    18956493 :                   t563 = rs**2
    1975    18956493 :                   t564 = 0.1e1_dp/t563
    1976    18956493 :                   t576 = t175**2
    1977    18956493 :                   t578 = t12/t576
    1978    18956493 :                   t579 = t27**2
    1979    18956493 :                   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    18956493 :                                     0.2e1_dp
    1987    18956493 :                   e_c_u_01rhoa = e_c_u_0rhoa
    1988    18956493 :                   t588 = alpha_1_2*rsrhoa
    1989    18956493 :                   t590 = t199*t211*t212
    1990    18956493 :                   t595 = t32/t198/t41
    1991    18956493 :                   t596 = t211**2
    1992    18956493 :                   t600 = beta_1_2*t548
    1993    18956493 :                   t606 = beta_3_2*t178
    1994    18956493 :                   t611 = t38**2
    1995    18956493 :                   t624 = t198**2
    1996    18956493 :                   t626 = t32/t624
    1997    18956493 :                   t627 = t45**2
    1998    18956493 :                   t628 = 0.1e1_dp/t627
    1999    18956493 :                   t636 = alpha_1_3*rsrhoa
    2000    18956493 :                   t638 = t220*t232*t233
    2001    18956493 :                   t643 = t50/t219/t59
    2002    18956493 :                   t644 = t232**2
    2003    18956493 :                   t648 = beta_1_3*t548
    2004    18956493 :                   t654 = beta_3_3*t178
    2005    18956493 :                   t659 = t56**2
    2006    18956493 :                   t672 = t219**2
    2007    18956493 :                   t674 = t50/t672
    2008    18956493 :                   t675 = t63**2
    2009    18956493 :                   t676 = 0.1e1_dp/t675
    2010    18956493 :                   alpha_c1rhoa = alpha_crhoa
    2011    18956493 :                   t681 = 0.1e1_dp/t87
    2012    18956493 :                   t682 = chirhoa**2
    2013    18956493 :                   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    18956493 :                                0.4e1_dp/0.3e1_dp*t74*chirhoarhoa)*t69
    2017    18956493 :                   f1rhoa = frhoa
    2018    18956493 :                   t705 = alpha_c1rhoa*f
    2019    18956493 :                   t708 = alpha_c*f1rhoa
    2020    18956493 :                   t711 = t78*t79
    2021    18956493 :                   t726 = e_c_u_1rhoa - e_c_u_01rhoa
    2022    18956493 :                   t733 = t726*f
    2023    18956493 :                   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    18956493 :                          + 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    18956493 :                                            + t745
    2046             :                   epsilon_c_unif1rhoa = e_c_u_01rhoa + t705*t82 + t708*t82 - &
    2047    18956493 :                                         t248 + t733*t80 + t736*t80 + t256
    2048    18956493 :                   t750 = 0.1e1_dp/t72
    2049    18956493 :                   t755 = 0.1e1_dp/t75
    2050             :                   phirhoarhoa = -t750*t682/0.9e1_dp + t257*chirhoarhoa/ &
    2051    18956493 :                                 0.3e1_dp - t755*t682/0.9e1_dp - t259*chirhoarhoa/0.3e1_dp
    2052    18956493 :                   phi1rhoa = phirhoa
    2053    18956493 :                   t763 = t90**2
    2054    18956493 :                   k_frhoarhoa = -0.2e1_dp/0.9e1_dp*t3/t262/t91*t763
    2055    18956493 :                   t767 = 0.1e1_dp/t94/t93
    2056    18956493 :                   t768 = k_frhoa**2
    2057    18956493 :                   k_s1rhoa = k_srhoa
    2058    18956493 :                   t775 = my_norm_drho*t105*t97
    2059    18956493 :                   t776 = t2*phirhoa
    2060    18956493 :                   t779 = t269*t273
    2061    18956493 :                   t785 = t269*t277*phirhoa/0.2e1_dp
    2062    18956493 :                   t789 = t2*k_srhoa
    2063    18956493 :                   t795 = t96/t272/k_s
    2064    18956493 :                   t798 = t273*t163
    2065    18956493 :                   t801 = t96*t798*k_srhoa/0.2e1_dp
    2066    18956493 :                   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    18956493 :                               /0.2e1_dp + t812
    2073             :                   t1rhoa = -t269*t98*phi1rhoa/0.2e1_dp - t96*t274*k_s1rhoa &
    2074    18956493 :                            /0.2e1_dp - t278/0.2e1_dp
    2075    18956493 :                   t820 = t101/t280/t108
    2076    18956493 :                   t821 = t107**2
    2077    18956493 :                   t822 = t289*t821
    2078    18956493 :                   t823 = epsilon_c_unif1rhoa*t100
    2079    18956493 :                   t825 = t285*phi1rhoa
    2080    18956493 :                   t828 = -t823*t105 + 0.3e1_dp*t102*t825
    2081    18956493 :                   t839 = 0.1e1_dp/t284/phi
    2082    18956493 :                   t840 = t839*phirhoa
    2083    18956493 :                   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    18956493 :                               t107
    2089    18956493 :                   A1rhoa = -t101*t281*t828*t107
    2090    18956493 :                   t858 = gamma_var*phi
    2091    18956493 :                   t865 = A1rhoa*t111
    2092    18956493 :                   t867 = 0.2e1_dp*t303*t1rhoa
    2093    18956493 :                   t868 = t865 + t867
    2094    18956493 :                   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    18956493 :                          - t310*t313*t876
    2097    18956493 :                   t880 = t879*t325
    2098    18956493 :                   t904 = t306*t119
    2099    18956493 :                   t908 = Arhoarhoa*t111
    2100    18956493 :                   t909 = Arhoa*t
    2101    18956493 :                   t911 = 0.2e1_dp*t909*t1rhoa
    2102    18956493 :                   t914 = 0.2e1_dp*A1rhoa*t*trhoa
    2103    18956493 :                   t917 = 0.2e1_dp*A*t1rhoa*trhoa
    2104    18956493 :                   t919 = 0.2e1_dp*t303*trhoarhoa
    2105    18956493 :                   t924 = t306*t312
    2106    18956493 :                   t936 = t113/t311/t118
    2107    18956493 :                   t944 = A*t317
    2108    18956493 :                   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    18956493 :                          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    18956493 :                          t321 + 0.2e1_dp*t310*t936*t321*t876 - t310*t313*t959
    2119    18956493 :                   t965 = t122**2
    2120    18956493 :                   t966 = 0.1e1_dp/t965
    2121    18956493 :                   t967 = t324*t966
    2122    18956493 :                   kf_arhoarhoa = -0.2e1_dp/0.9e1_dp*t124/t329/t125*t763
    2123    18956493 :                   ex_unif_a1rhoa = ex_unif_arhoa
    2124    18956493 :                   t985 = kf_arhoa**2
    2125    18956493 :                   s_a1rhoa = s_arhoa
    2126    18956493 :                   t998 = mu**2
    2127    18956493 :                   t999 = 0.1e1_dp/t344/t137*t998
    2128    18956493 :                   t1000 = t999*t133
    2129    18956493 :                   t1001 = s_arhoa*t135
    2130    18956493 :                   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    18956493 :                                                                                                                   - t110*t967*t879))
    2147             : 
    2148    18956493 :                   chirhoarhob = 0.2e1_dp*t519
    2149    18956493 :                   rsrhoarhob = rsrhoarhoa
    2150    18956493 :                   t1031 = t176*t368*t191
    2151    18956493 :                   t1033 = alpha_1_1*rsrhob
    2152    18956493 :                   t1038 = rsrhoa*rsrhob
    2153    18956493 :                   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    18956493 :                                     t14*t368/0.2e1_dp
    2161    18956493 :                   t1069 = t199*t382*t212
    2162    18956493 :                   t1071 = alpha_1_2*rsrhob
    2163    18956493 :                   t1104 = t220*t396*t233
    2164    18956493 :                   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    18956493 :                               t69
    2169    18956493 :                   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    18956493 :                           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    18956493 :                                            t1193
    2192             :                   phirhoarhob = -t750*chirhoa*chirhob/0.9e1_dp + t257* &
    2193             :                                 chirhoarhob/0.3e1_dp - t755*chirhoa*chirhob/0.9e1_dp - t259 &
    2194    18956493 :                                 *chirhoarhob/0.3e1_dp
    2195    18956493 :                   k_frhoarhob = k_frhoarhoa
    2196    18956493 :                   t1228 = t269*t277*phirhob/0.2e1_dp
    2197    18956493 :                   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    18956493 :                               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    18956493 :                               t107
    2208    18956493 :                   t1269 = t445*t119
    2209    18956493 :                   t1283 = Arhoarhob*t111
    2210    18956493 :                   t1285 = 0.2e1_dp*t909*trhob
    2211    18956493 :                   t1286 = Arhob*t
    2212    18956493 :                   t1288 = 0.2e1_dp*t1286*trhoa
    2213    18956493 :                   t1291 = 0.2e1_dp*A*trhob*trhoa
    2214    18956493 :                   t1293 = 0.2e1_dp*t303*trhoarhob
    2215    18956493 :                   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    18956493 :                           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    18956493 :                           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    18956493 :                                                   - t110*t967*t456))
    2233             : 
    2234    18956493 :                   chirhobrhob = 0.2e1_dp*t163 + 0.2e1_dp*t519
    2235    18956493 :                   rsrhobrhob = rsrhoarhob
    2236    18956493 :                   t1342 = t368**2
    2237    18956493 :                   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    18956493 :                                     t1342*t580*t14/0.2e1_dp
    2245    18956493 :                   e_c_u_01rhob = e_c_u_0rhob
    2246    18956493 :                   t1377 = t382**2
    2247    18956493 :                   t1411 = t396**2
    2248    18956493 :                   alpha_c1rhob = alpha_crhob
    2249    18956493 :                   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    18956493 :                                0.4e1_dp/0.3e1_dp*t74*chirhobrhob)*t69
    2253    18956493 :                   f1rhob = frhob
    2254    18956493 :                   t1462 = alpha_c1rhob*f
    2255    18956493 :                   t1465 = alpha_c*f1rhob
    2256    18956493 :                   t1482 = e_c_u_1rhob - e_c_u_01rhob
    2257    18956493 :                   t1489 = t1482*f
    2258    18956493 :                   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    18956493 :                           *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    18956493 :                                            *t711*t1440 + t1501
    2281             :                   epsilon_c_unif1rhob = e_c_u_01rhob + t1462*t82 + t1465*t82 - &
    2282    18956493 :                                         t409 + t1489*t80 + t1492*t80 + t417
    2283             :                   phirhobrhob = -t750*t1440/0.9e1_dp + t257*chirhobrhob/ &
    2284    18956493 :                                 0.3e1_dp - t755*t1440/0.9e1_dp - t259*chirhobrhob/0.3e1_dp
    2285    18956493 :                   phi1rhob = phirhob
    2286    18956493 :                   t1514 = k_frhob**2
    2287    18956493 :                   k_s1rhob = k_srhob
    2288    18956493 :                   t1520 = t2*phirhob
    2289    18956493 :                   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    18956493 :                               k_s1rhob/0.2e1_dp + t812
    2296             :                   t1rhob = -t269*t98*phi1rhob/0.2e1_dp - t96*t274*k_s1rhob &
    2297    18956493 :                            /0.2e1_dp - t278/0.2e1_dp
    2298    18956493 :                   t1550 = epsilon_c_unif1rhob*t100
    2299    18956493 :                   t1552 = t285*phi1rhob
    2300    18956493 :                   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    18956493 :                               t432*t1555*t107
    2306    18956493 :                   A1rhob = -t101*t281*t1555*t107
    2307    18956493 :                   t1588 = A1rhob*t111
    2308    18956493 :                   t1590 = 0.2e1_dp*t303*t1rhob
    2309    18956493 :                   t1591 = t1588 + t1590
    2310             :                   t1599 = t1588 + t1590 + 0.2e1_dp*t314*A1rhob + 0.4e1_dp*t318 &
    2311    18956493 :                           *t1rhob
    2312             :                   t1602 = 0.2e1_dp*t297*t298*t1rhob + t101*t111*t1591* &
    2313    18956493 :                           t119 - t310*t313*t1599
    2314    18956493 :                   t1603 = t1602*t325
    2315    18956493 :                   t1630 = Arhobrhob*t111
    2316    18956493 :                   t1632 = 0.2e1_dp*t1286*t1rhob
    2317    18956493 :                   t1635 = 0.2e1_dp*A1rhob*t*trhob
    2318    18956493 :                   t1638 = 0.2e1_dp*A*t1rhob*trhob
    2319    18956493 :                   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    18956493 :                           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    18956493 :                           t313*t1674
    2331    18956493 :                   t1680 = t456*t966
    2332    18956493 :                   kf_brhobrhob = -0.2e1_dp/0.9e1_dp*t124/t460/t142*t763
    2333    18956493 :                   ex_unif_b1rhob = ex_unif_brhob
    2334    18956493 :                   t1698 = kf_brhob**2
    2335    18956493 :                   s_b1rhob = s_brhob
    2336    18956493 :                   t1711 = 0.1e1_dp/t475/t153*t998
    2337    18956493 :                   t1712 = t1711*t150
    2338    18956493 :                   t1713 = s_brhob*t135
    2339    18956493 :                   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    18956493 :                                                                                                            t325 - t110*t1680*t1602))
    2356    18956493 :                   t1739 = t268*t97
    2357    18956493 :                   t1741 = t95*t273
    2358    18956493 :                   t1743 = t488*t163
    2359             :                   trhoanorm_drho = -t1739*t776/0.2e1_dp - t1741*t789/ &
    2360    18956493 :                                    0.2e1_dp - t1743/0.2e1_dp
    2361    18956493 :                   t1748 = t101*tnorm_drho
    2362    18956493 :                   t1765 = t909*tnorm_drho
    2363    18956493 :                   t1766 = t494*trhoa
    2364    18956493 :                   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    18956493 :                                                        trhoanorm_drho)
    2375             : 
    2376             :                   e_ra_ndr(ii) = e_ra_ndr(ii) + &
    2377             :                                  scale_ec*(Hnorm_drho + my_rho*(0.3e1_dp* &
    2378    18956493 :                                                                 t293*t506*phirhoa + t110*t1801*t325 - t110*t967*t505))
    2379             : 
    2380             :                   trhobnorm_drho = -t1739*t1520/0.2e1_dp - t1741*t1529/ &
    2381    18956493 :                                    0.2e1_dp - t1743/0.2e1_dp
    2382    18956493 :                   t1829 = t1286*tnorm_drho
    2383    18956493 :                   t1830 = t494*trhob
    2384    18956493 :                   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    18956493 :                                                        trhobnorm_drho)
    2395             : 
    2396             :                   e_rb_ndr(ii) = e_rb_ndr(ii) + &
    2397             :                                  scale_ec*(Hnorm_drho + my_rho*(0.3e1_dp* &
    2398    18956493 :                                                                 t293*t506*phirhob + t110*t1865*t325 - t110*t1680*t505))
    2399             : 
    2400    18956493 :                   t1871 = tnorm_drho**2
    2401    18956493 :                   t1876 = A*t1871
    2402    18956493 :                   t1888 = t502**2
    2403    18956493 :                   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    18956493 :                                                    *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    18956493 :                                                                                   kf_arhoa/0.2e1_dp - t129*t341/0.2e1_dp)))/0.2e1_dp
    2419             : 
    2420    18956493 :                   t1922 = s_anorm_drhoa**2
    2421             : 
    2422             :                   e_ndra_ndra(ii) = e_ndra_ndra(ii) + &
    2423             :                                     scale_ex*t140*(-0.8e1_dp*t999* &
    2424    18956493 :                                                    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    18956493 :                                                                                   kf_brhob/0.2e1_dp - t146*t472/0.2e1_dp)))/0.2e1_dp
    2432             : 
    2433    18956493 :                   t1949 = s_bnorm_drhob**2
    2434             :                   e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii) + &
    2435             :                                     scale_ex*t156*(-0.8e1_dp*t1711* &
    2436    18956493 :                                                    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       18025 :    END SUBROUTINE pbe_lsd_calc
    2443             : 
    2444             : END MODULE xc_pbe

Generated by: LCOV version 1.15