LCOV - code coverage report
Current view: top level - src/xc - xc_ke_gga.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 48.0 % 252 121
Test Date: 2025-07-25 12:55:17 Functions: 73.3 % 15 11

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Calculate the several different kinetic energy functionals
      10              : !>      with a GGA form
      11              : !> \par History
      12              : !>      JGH (26.02.2003) : OpenMP enabled
      13              : !>      fawzi (04.2004)  : adapted to the new xc interface
      14              : !> \author JGH (20.02.2002)
      15              : ! **************************************************************************************************
      16              : MODULE xc_ke_gga
      17              : 
      18              :    USE cp_array_utils,                  ONLY: cp_3d_r_cp_type
      19              :    USE kinds,                           ONLY: dp
      20              :    USE mathconstants,                   ONLY: pi
      21              :    USE xc_derivative_desc,              ONLY: deriv_norm_drho,&
      22              :                                               deriv_norm_drhoa,&
      23              :                                               deriv_norm_drhob,&
      24              :                                               deriv_rho,&
      25              :                                               deriv_rhoa,&
      26              :                                               deriv_rhob
      27              :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      28              :                                               xc_dset_get_derivative
      29              :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      30              :                                               xc_derivative_type
      31              :    USE xc_functionals_utilities,        ONLY: calc_wave_vector,&
      32              :                                               set_util
      33              :    USE xc_input_constants,              ONLY: ke_lc,&
      34              :                                               ke_llp,&
      35              :                                               ke_ol1,&
      36              :                                               ke_ol2,&
      37              :                                               ke_pbe,&
      38              :                                               ke_pw86,&
      39              :                                               ke_pw91,&
      40              :                                               ke_t92
      41              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      42              :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      43              :                                               xc_rho_set_type
      44              : #include "../base/base_uses.f90"
      45              : 
      46              :    IMPLICIT NONE
      47              : 
      48              :    PRIVATE
      49              : 
      50              :    REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
      51              :                                f23 = 2.0_dp*f13, &
      52              :                                f43 = 4.0_dp*f13, &
      53              :                                f53 = 5.0_dp*f13
      54              : 
      55              :    PUBLIC :: ke_gga_info, ke_gga_lda_eval, ke_gga_lsd_eval
      56              : 
      57              :    REAL(KIND=dp) :: cf, b, b_lda, b_lsd, flda, flsd, sfac, t13
      58              :    REAL(KIND=dp) :: fact, tact
      59              :    REAL(KIND=dp) :: eps_rho
      60              : 
      61              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_ke_gga'
      62              : 
      63              : CONTAINS
      64              : 
      65              : ! **************************************************************************************************
      66              : !> \brief ...
      67              : !> \param cutoff ...
      68              : ! **************************************************************************************************
      69         1648 :    SUBROUTINE ke_gga_init(cutoff)
      70              : 
      71              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
      72              : 
      73         1648 :       eps_rho = cutoff
      74         1648 :       CALL set_util(cutoff)
      75              : 
      76         1648 :       cf = 0.3_dp*(3.0_dp*pi*pi)**f23
      77         1648 :       flda = cf
      78         1648 :       flsd = flda*2.0_dp**f23
      79              : !   the_factor 2^(1/3) for LDA is here
      80         1648 :       b_lda = 2.0_dp**f43*(3.0_dp*pi*pi)**(f13)
      81         1648 :       b_lsd = 2.0_dp*(3.0_dp*pi*pi)**(f13)
      82         1648 :       sfac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
      83         1648 :       t13 = 2.0_dp**f13
      84              : 
      85         1648 :    END SUBROUTINE ke_gga_init
      86              : 
      87              : ! **************************************************************************************************
      88              : !> \brief ...
      89              : !> \param functional ...
      90              : !> \param lsd ...
      91              : !> \param reference ...
      92              : !> \param shortform ...
      93              : !> \param needs ...
      94              : !> \param max_deriv ...
      95              : ! **************************************************************************************************
      96         1786 :    SUBROUTINE ke_gga_info(functional, lsd, reference, shortform, needs, max_deriv)
      97              :       INTEGER, INTENT(in)                                :: functional
      98              :       LOGICAL, INTENT(in)                                :: lsd
      99              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
     100              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
     101              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
     102              : 
     103         1786 :       IF (PRESENT(reference)) THEN
     104            0 :          SELECT CASE (functional)
     105              :          CASE (ke_ol1)
     106              :             reference = "H. Ou-Yang and M. Levy, "// &
     107            0 :                         "Intl. J. Quant. Chem. 40, 379 (1991); Functional 1"
     108              :          CASE (ke_ol2)
     109              :             reference = "H. Ou-Yang and M. Levy, "// &
     110            0 :                         "Intl. J. Quant. Chem. 40, 379 (1991); Functional 2"
     111              :          CASE (ke_llp)
     112            0 :             reference = "H. Lee, C. Lee, R.G. Parr, Phys. Rev. A, 44, 768 (1991)"
     113              :          CASE (ke_pw86)
     114            0 :             reference = "J.P. Perdew and Y. Wang, Phys. Rev. B, 33, 8800 (1986)"
     115              :          CASE (ke_pw91)
     116            0 :             reference = "J.P. Perdew and Y. Wang, Electronic Structure of Solids 91"
     117              :          CASE (ke_lc)
     118            0 :             reference = "A. Lembarki and H. Chermette, Phys. Rev. A, 50, 5328 (1994)"
     119              :          CASE (ke_t92)
     120            0 :             reference = "A.J. Thakkar, Phys. Rev. A, 46, 6920 (1992)"
     121              :          CASE (ke_pbe)
     122            0 :             reference = "J.P.Perdew, K.Burke, M.Ernzerhof, Phys. Rev. Letter, 77, 3865 (1996)"
     123              :          END SELECT
     124            0 :          IF (.NOT. lsd) THEN
     125            0 :             IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
     126            0 :                reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {spin unpolarized}'
     127              :             END IF
     128              :          END IF
     129              :       END IF
     130         1786 :       IF (PRESENT(shortform)) THEN
     131            0 :          SELECT CASE (functional)
     132              :          CASE (ke_ol1)
     133            0 :             shortform = "Ou-Yang-Levy Functional 1"
     134              :          CASE (ke_ol2)
     135            0 :             shortform = "Ou-Yang-Levy Functional 2"
     136              :          CASE (ke_llp)
     137            0 :             shortform = "Lee-Lee-Parr Functional"
     138              :          CASE (ke_pw86)
     139            0 :             shortform = "Perdew-Wang 1986 Functional (kinetic energy)"
     140              :          CASE (ke_pw91)
     141            0 :             shortform = "Perdew-Wang 1991 Functional (kinetic energy)"
     142              :          CASE (ke_lc)
     143            0 :             shortform = "Lembarki-Chermette kinetic energy functional"
     144              :          CASE (ke_t92)
     145            0 :             shortform = "Thakkar 1992 Functional"
     146              :          CASE (ke_pbe)
     147            0 :             shortform = "Perdew-Burke-Ernzerhof Functional (kinetic energy)"
     148              :          END SELECT
     149            0 :          IF (.NOT. lsd) THEN
     150            0 :             IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
     151            0 :                shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {spin unpolarized}'
     152              :             END IF
     153              :          END IF
     154              :       END IF
     155         1786 :       IF (PRESENT(needs)) THEN
     156         1786 :          IF (lsd) THEN
     157            0 :             needs%rho_spin = .TRUE.
     158            0 :             needs%rho_spin_1_3 = .TRUE.
     159            0 :             needs%norm_drho_spin = .TRUE.
     160              :          ELSE
     161         1786 :             needs%rho = .TRUE.
     162         1786 :             needs%rho_1_3 = .TRUE.
     163         1786 :             needs%norm_drho = .TRUE.
     164              :          END IF
     165              :       END IF
     166         1786 :       IF (PRESENT(max_deriv)) max_deriv = 3
     167              : 
     168         1786 :    END SUBROUTINE ke_gga_info
     169              : 
     170              : ! **************************************************************************************************
     171              : !> \brief ...
     172              : !> \param functional ...
     173              : !> \param rho_set ...
     174              : !> \param deriv_set ...
     175              : !> \param order ...
     176              : ! **************************************************************************************************
     177         1648 :    SUBROUTINE ke_gga_lda_eval(functional, rho_set, deriv_set, order)
     178              : 
     179              :       INTEGER, INTENT(IN)                                :: functional
     180              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     181              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     182              :       INTEGER, INTENT(IN)                                :: order
     183              : 
     184              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'ke_gga_lda_eval'
     185              : 
     186              :       INTEGER                                            :: handle, m, npoints
     187              :       INTEGER, DIMENSION(2, 3)                           :: bo
     188              :       REAL(KIND=dp)                                      :: drho_cutoff, rho_cutoff
     189              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: s
     190              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: fs
     191         1648 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
     192         1648 :          e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
     193         1648 :          e_rho_rho_rho, grho, rho, rho13
     194              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     195              : 
     196         1648 :       CALL timeset(routineN, handle)
     197         1648 :       NULLIFY (rho, rho13, e_0, e_rho, e_ndrho, &
     198         1648 :                e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
     199         1648 :                e_rho_rho_rho, e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho)
     200         1648 :       m = ABS(order)
     201              : 
     202              :       CALL xc_rho_set_get(rho_set, rho_1_3=rho13, rho=rho, &
     203              :                           norm_drho=grho, local_bounds=bo, rho_cutoff=rho_cutoff, &
     204         1648 :                           drho_cutoff=drho_cutoff)
     205         1648 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     206         1648 :       CALL ke_gga_init(rho_cutoff)
     207              : 
     208         4944 :       ALLOCATE (s(npoints))
     209         6592 :       ALLOCATE (fs(npoints, m + 1))
     210              : 
     211              : !      s = norm_drho/(rho^(4/3)*2*(pi*pi*3)^(1/3))
     212         1648 :       CALL calc_wave_vector("p", rho, grho, s)
     213              : 
     214         1648 :       fact = flda
     215              : !      Definition of s has changed
     216         1648 :       b = b_lda
     217              : !       tact = t13
     218         1648 :       tact = 1.0_dp
     219              : 
     220         1648 :       SELECT CASE (functional)
     221              :       CASE (ke_ol1)
     222            0 :          CALL efactor_ol1(s, fs, m)
     223            0 :          CPABORT("OL1 functional currently not working properly")
     224              :       CASE (ke_ol2)
     225            0 :          CALL efactor_ol2(s, fs, m)
     226            0 :          CPABORT("OL2 functional currently not working properly")
     227              :       CASE (ke_llp)
     228          736 :          CALL efactor_llp(s, fs, m)
     229              :       CASE (ke_pw86)
     230           54 :          CALL efactor_pw86(s, fs, m)
     231              :       CASE (ke_pw91)
     232           54 :          CALL efactor_pw91(s, fs, m, 1)
     233              :       CASE (ke_lc)
     234           54 :          CALL efactor_pw91(s, fs, m, 2)
     235              :       CASE (ke_t92)
     236          174 :          CALL efactor_t92(s, fs, m)
     237              :       CASE (ke_pbe)
     238          576 :          CALL efactor_pbex(s, fs, m, 1)
     239              :       CASE DEFAULT
     240         1648 :          CPABORT("Unsupported functional")
     241              :       END SELECT
     242              : 
     243         1648 :       IF (order >= 0) THEN
     244              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     245         1648 :                                          allocate_deriv=.TRUE.)
     246         1648 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     247              : 
     248         1648 :          CALL kex_p_0(rho, rho13, fs, e_0, npoints)
     249              :       END IF
     250              : 
     251         1648 :       IF (order >= 1 .OR. order == -1) THEN
     252              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     253         1648 :                                          allocate_deriv=.TRUE.)
     254         1648 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     255              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     256         1648 :                                          allocate_deriv=.TRUE.)
     257         1648 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     258              : 
     259              :          CALL kex_p_1(rho, rho13, s, fs, e_rho=e_rho, e_ndrho=e_ndrho, &
     260         1648 :                       npoints=npoints)
     261              :       END IF
     262         1648 :       IF (order >= 2 .OR. order == -2) THEN
     263              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     264          138 :                                          allocate_deriv=.TRUE.)
     265          138 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     266              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
     267          138 :                                          allocate_deriv=.TRUE.)
     268          138 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
     269              :          deriv => xc_dset_get_derivative(deriv_set, &
     270          138 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     271          138 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     272              : 
     273              :          CALL kex_p_2(rho, rho13, s, fs, e_rho_rho=e_rho_rho, &
     274          138 :                       e_rho_ndrho=e_rho_ndrho, e_ndrho_ndrho=e_ndrho_ndrho, npoints=npoints)
     275              :       END IF
     276         1648 :       IF (order >= 3 .OR. order == -3) THEN
     277              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
     278            0 :                                          allocate_deriv=.TRUE.)
     279            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     280              :          deriv => xc_dset_get_derivative(deriv_set, &
     281            0 :                                          [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.TRUE.)
     282            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
     283              :          deriv => xc_dset_get_derivative(deriv_set, &
     284            0 :                                          [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     285            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
     286              :          deriv => xc_dset_get_derivative(deriv_set, &
     287            0 :                                          [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     288            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
     289              : 
     290              :          CALL kex_p_3(rho, rho13, s, fs, e_rho_rho_rho=e_rho_rho_rho, &
     291              :                       e_rho_rho_ndrho=e_rho_rho_ndrho, e_rho_ndrho_ndrho=e_rho_ndrho_ndrho, &
     292            0 :                       e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, npoints=npoints)
     293              :       END IF
     294         1648 :       IF (order > 3 .OR. order < -3) THEN
     295            0 :          CPABORT("derivatives bigger than 3 not implemented")
     296              :       END IF
     297              : 
     298         1648 :       DEALLOCATE (s)
     299         1648 :       DEALLOCATE (fs)
     300              : 
     301         1648 :       CALL timestop(handle)
     302              : 
     303         1648 :    END SUBROUTINE ke_gga_lda_eval
     304              : 
     305              : ! **************************************************************************************************
     306              : !> \brief ...
     307              : !> \param functional ...
     308              : !> \param rho_set ...
     309              : !> \param deriv_set ...
     310              : !> \param order ...
     311              : ! **************************************************************************************************
     312            0 :    SUBROUTINE ke_gga_lsd_eval(functional, rho_set, deriv_set, order)
     313              : 
     314              :       INTEGER, INTENT(IN)                                :: functional
     315              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     316              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     317              :       INTEGER, INTENT(IN), OPTIONAL                      :: order
     318              : 
     319              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ke_gga_lsd_eval'
     320              :       INTEGER, DIMENSION(2), PARAMETER :: &
     321              :          norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob], &
     322              :          rho_spin_name = [deriv_rhoa, deriv_rhob]
     323              : 
     324              :       INTEGER                                            :: handle, ispin, m, npoints
     325              :       INTEGER, DIMENSION(2, 3)                           :: bo
     326              :       REAL(KIND=dp)                                      :: rho_cutoff
     327            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: s
     328            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: fs
     329            0 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
     330            0 :          e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
     331            0 :          e_rho_rho_rho
     332            0 :       TYPE(cp_3d_r_cp_type), DIMENSION(2)                :: norm_drho, rho, rho_1_3
     333              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     334              : 
     335            0 :       CALL timeset(routineN, handle)
     336            0 :       NULLIFY (e_0, e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_rho_ndrho_ndrho, &
     337            0 :                e_rho_ndrho, e_rho_rho_ndrho, e_rho, e_rho_rho, e_rho_rho_rho)
     338            0 :       DO ispin = 1, 2
     339            0 :          NULLIFY (norm_drho(ispin)%array, rho(ispin)%array, rho_1_3(ispin)%array)
     340              :       END DO
     341              : 
     342              :       CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
     343              :                           rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
     344              :                           rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
     345              :                           norm_drhob=norm_drho(2)%array, rho_cutoff=rho_cutoff, &
     346            0 :                           local_bounds=bo)
     347            0 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     348            0 :       m = ABS(order)
     349            0 :       CALL ke_gga_init(rho_cutoff)
     350              : 
     351            0 :       ALLOCATE (s(npoints))
     352            0 :       ALLOCATE (fs(npoints, m + 1))
     353              : 
     354            0 :       fact = flsd
     355            0 :       b = b_lsd
     356            0 :       tact = 1.0_dp
     357              : 
     358            0 :       DO ispin = 1, 2
     359              : 
     360            0 :          CALL calc_wave_vector("p", rho(ispin)%array, norm_drho(ispin)%array, s)
     361              : 
     362            0 :          SELECT CASE (functional)
     363              :          CASE (ke_ol1)
     364            0 :             CALL efactor_ol1(s, fs, m)
     365              :          CASE (ke_ol2)
     366            0 :             CALL efactor_ol2(s, fs, m)
     367              :          CASE (ke_llp)
     368            0 :             CALL efactor_llp(s, fs, m)
     369              :          CASE (ke_pw86)
     370            0 :             tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
     371            0 :             CALL efactor_pw86(s, fs, m, f2_lsd=tact)
     372            0 :             tact = 1.0_dp
     373              :          CASE (ke_pbe)
     374            0 :             tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
     375            0 :             CALL efactor_pbex(s, fs, m, 1, f2_lsd=tact)
     376            0 :             tact = 1.0_dp
     377              :          CASE (ke_pw91)
     378            0 :             tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
     379            0 :             CALL efactor_pw91(s, fs, m, 1, f2_lsd=tact)
     380            0 :             tact = 1.0_dp
     381              :          CASE (ke_lc)
     382            0 :             tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
     383            0 :             CALL efactor_pw91(s, fs, m, 2, f2_lsd=tact)
     384            0 :             tact = 1.0_dp
     385              :          CASE (ke_t92)
     386            0 :             CALL efactor_t92(s, fs, m)
     387              :          CASE DEFAULT
     388            0 :             CPABORT("Unsupported functional")
     389              :          END SELECT
     390              : 
     391            0 :          IF (order >= 0) THEN
     392              :             deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     393            0 :                                             allocate_deriv=.TRUE.)
     394            0 :             CALL xc_derivative_get(deriv, deriv_data=e_0)
     395              : 
     396              :             CALL kex_p_0(rho(ispin)%array, rho_1_3(ispin)%array, fs, &
     397            0 :                          e_0=e_0, npoints=npoints)
     398              :          END IF
     399            0 :          IF (order >= 1 .OR. order == -1) THEN
     400              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
     401            0 :                                             allocate_deriv=.TRUE.)
     402            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho)
     403              :             deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
     404            0 :                                             allocate_deriv=.TRUE.)
     405            0 :             CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     406              : 
     407              :             CALL kex_p_1(rho=rho(ispin)%array, &
     408              :                          r13=rho_1_3(ispin)%array, s=s, fs=fs, e_rho=e_rho, &
     409            0 :                          e_ndrho=e_ndrho, npoints=npoints)
     410              :          END IF
     411            0 :          IF (order >= 2 .OR. order == -2) THEN
     412              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     413            0 :                                                         rho_spin_name(ispin)], allocate_deriv=.TRUE.)
     414            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     415              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     416            0 :                                                         norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
     417            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
     418              :             deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
     419            0 :                                                         norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
     420            0 :             CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     421              : 
     422              :             CALL kex_p_2(rho(ispin)%array, rho_1_3(ispin)%array, &
     423            0 :                          s, fs, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, npoints)
     424              :          END IF
     425            0 :          IF (order >= 3 .OR. order == -3) THEN
     426              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     427              :                                                         rho_spin_name(ispin), rho_spin_name(ispin)], &
     428            0 :                                             allocate_deriv=.TRUE.)
     429            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     430              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     431              :                                                         rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
     432            0 :                                             allocate_deriv=.TRUE.)
     433            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
     434              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     435              :                                                         norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
     436            0 :                                             allocate_deriv=.TRUE.)
     437            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
     438              :             deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
     439              :                                                         norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
     440            0 :                                             allocate_deriv=.TRUE.)
     441            0 :             CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
     442              : 
     443              :             CALL kex_p_3(rho(ispin)%array, &
     444              :                          rho_1_3(ispin)%array, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
     445            0 :                          e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
     446              :          END IF
     447            0 :          IF (order > 3 .OR. order < -3) THEN
     448            0 :             CPABORT("derivatives bigger than 3 not implemented")
     449              :          END IF
     450              : 
     451              :       END DO
     452              : 
     453            0 :       DEALLOCATE (s)
     454            0 :       DEALLOCATE (fs)
     455            0 :       CALL timestop(handle)
     456            0 :    END SUBROUTINE ke_gga_lsd_eval
     457              : 
     458              : ! **************************************************************************************************
     459              : !> \brief ...
     460              : !> \param rho ...
     461              : !> \param r13 ...
     462              : !> \param fs ...
     463              : !> \param e_0 ...
     464              : !> \param npoints ...
     465              : ! **************************************************************************************************
     466         1648 :    SUBROUTINE kex_p_0(rho, r13, fs, e_0, npoints)
     467              : 
     468              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13
     469              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: fs
     470              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_0
     471              :       INTEGER, INTENT(in)                                :: npoints
     472              : 
     473              :       INTEGER                                            :: ip
     474              : 
     475              : !$OMP     PARALLEL DO DEFAULT(NONE) &
     476              : !$OMP                 SHARED(npoints, rho, e_0, fact, r13, fs, eps_rho) &
     477         1648 : !$OMP                 PRIVATE(ip)
     478              : 
     479              :       DO ip = 1, npoints
     480              : 
     481              :          IF (rho(ip) > eps_rho) THEN
     482              :             e_0(ip) = e_0(ip) + fact*r13(ip)*r13(ip)*rho(ip)*fs(ip, 1)
     483              :          END IF
     484              : 
     485              :       END DO
     486              : 
     487              : !$OMP     END PARALLEL DO
     488              : 
     489         1648 :    END SUBROUTINE kex_p_0
     490              : 
     491              : ! **************************************************************************************************
     492              : !> \brief ...
     493              : !> \param rho ...
     494              : !> \param r13 ...
     495              : !> \param s ...
     496              : !> \param fs ...
     497              : !> \param e_rho ...
     498              : !> \param e_ndrho ...
     499              : !> \param npoints ...
     500              : ! **************************************************************************************************
     501         1648 :    SUBROUTINE kex_p_1(rho, r13, s, fs, e_rho, e_ndrho, npoints)
     502              : 
     503              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13, s
     504              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: fs
     505              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho, e_ndrho
     506              :       INTEGER, INTENT(in)                                :: npoints
     507              : 
     508              :       INTEGER                                            :: ip
     509              :       REAL(KIND=dp)                                      :: a0, a1, sx, sy
     510              : 
     511              : !$OMP     PARALLEL DO DEFAULT(NONE) &
     512              : !$OMP                 SHARED(npoints, rho, eps_rho, fact, r13, sfac, tact) &
     513              : !$OMP                 SHARED(fs, e_rho, e_ndrho, s) &
     514         1648 : !$OMP                 PRIVATE(ip,a0,a1,sx,sy)
     515              : 
     516              :       DO ip = 1, npoints
     517              : 
     518              :          IF (rho(ip) > eps_rho) THEN
     519              : 
     520              :             a0 = fact*r13(ip)*r13(ip)*rho(ip)
     521              :             a1 = f53*fact*r13(ip)*r13(ip)
     522              :             sx = -f43*s(ip)/rho(ip)
     523              :             sy = sfac*tact/(r13(ip)*rho(ip))
     524              :             e_rho(ip) = e_rho(ip) + a1*fs(ip, 1) + a0*fs(ip, 2)*sx
     525              :             e_ndrho(ip) = e_ndrho(ip) + a0*fs(ip, 2)*sy
     526              : 
     527              :          END IF
     528              : 
     529              :       END DO
     530              : 
     531              : !$OMP     END PARALLEL DO
     532              : 
     533         1648 :    END SUBROUTINE kex_p_1
     534              : 
     535              : ! **************************************************************************************************
     536              : !> \brief ...
     537              : !> \param rho ...
     538              : !> \param r13 ...
     539              : !> \param s ...
     540              : !> \param fs ...
     541              : !> \param e_rho_rho ...
     542              : !> \param e_rho_ndrho ...
     543              : !> \param e_ndrho_ndrho ...
     544              : !> \param npoints ...
     545              : ! **************************************************************************************************
     546          138 :    SUBROUTINE kex_p_2(rho, r13, s, fs, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
     547              :                       npoints)
     548              : 
     549              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13, s
     550              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: fs
     551              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
     552              :       INTEGER, INTENT(in)                                :: npoints
     553              : 
     554              :       INTEGER                                            :: ip
     555              :       REAL(KIND=dp)                                      :: a0, a1, a2, sx, sxx, sxy, sy
     556              : 
     557              : !$OMP     PARALLEL DO DEFAULT(NONE) &
     558              : !$OMP                 SHARED (npoints, rho, eps_rho, fact, r13) &
     559              : !$OMP                 SHARED (e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, fs) &
     560              : !$OMP                 SHARED(s, sfac, tact) &
     561          138 : !$OMP                 PRIVATE(ip,a0,a1,a2,sx,sy,sxx,sxy)
     562              : 
     563              :       DO ip = 1, npoints
     564              : 
     565              :          IF (rho(ip) > eps_rho) THEN
     566              : 
     567              :             a0 = fact*r13(ip)*r13(ip)*rho(ip)
     568              :             a1 = f53*fact*r13(ip)*r13(ip)
     569              :             a2 = f23*f53*fact/r13(ip)
     570              :             sx = -f43*s(ip)/rho(ip)
     571              :             sy = sfac*tact/(r13(ip)*rho(ip))
     572              :             sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
     573              :             sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
     574              :             e_rho_rho(ip) = e_rho_rho(ip) + a2*fs(ip, 1) + 2.0_dp*a1*fs(ip, 2)*sx + &
     575              :                             a0*fs(ip, 3)*sx*sx + a0*fs(ip, 2)*sxx
     576              :             e_rho_ndrho(ip) = e_rho_ndrho(ip) + a1*fs(ip, 2)*sy + a0*fs(ip, 3)*sx*sy + &
     577              :                               a0*fs(ip, 2)*sxy
     578              :             e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + a0*fs(ip, 3)*sy*sy
     579              : 
     580              :          END IF
     581              : 
     582              :       END DO
     583              : 
     584              : !$OMP     END PARALLEL DO
     585              : 
     586          138 :    END SUBROUTINE kex_p_2
     587              : 
     588              : ! **************************************************************************************************
     589              : !> \brief ...
     590              : !> \param rho ...
     591              : !> \param r13 ...
     592              : !> \param s ...
     593              : !> \param fs ...
     594              : !> \param e_rho_rho_rho ...
     595              : !> \param e_rho_rho_ndrho ...
     596              : !> \param e_rho_ndrho_ndrho ...
     597              : !> \param e_ndrho_ndrho_ndrho ...
     598              : !> \param npoints ...
     599              : ! **************************************************************************************************
     600            0 :    SUBROUTINE kex_p_3(rho, r13, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
     601              :                       e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
     602              : 
     603              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13, s
     604              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: fs
     605              :       REAL(KIND=dp), DIMENSION(*), INTENT(inout)         :: e_rho_rho_rho, e_rho_rho_ndrho, &
     606              :                                                             e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
     607              :       INTEGER, INTENT(in)                                :: npoints
     608              : 
     609              :       INTEGER                                            :: ip
     610              :       REAL(KIND=dp)                                      :: a0, a1, a2, a3, sx, sxx, sxxx, sxxy, &
     611              :                                                             sxy, sy
     612              : 
     613              : !$OMP     PARALLEL DO DEFAULT(NONE) &
     614              : !$OMP                 SHARED(npoints, rho, eps_rho, fact, r13) &
     615              : !$OMP                 SHARED(s, sfac, tact, fs, e_rho_rho_rho) &
     616              : !$OMP                 SHARED(e_rho_rho_ndrho, e_rho_ndrho_ndrho) &
     617              : !$OMP                 SHARED(e_ndrho_ndrho_ndrho) &
     618            0 : !$OMP                 PRIVATE(ip,a0,a1,a2,a3,sx,sy,sxx,sxy,sxxx,sxxy)
     619              : 
     620              :       DO ip = 1, npoints
     621              : 
     622              :          IF (rho(ip) > eps_rho) THEN
     623              : 
     624              :             a0 = fact*r13(ip)*r13(ip)*rho(ip)
     625              :             a1 = f53*fact*r13(ip)*r13(ip)
     626              :             a2 = f23*f53*fact/r13(ip)
     627              :             a3 = -f13*f23*f53*fact/(r13(ip)*rho(ip))
     628              :             sx = -f43*s(ip)/rho(ip)
     629              :             sy = sfac*tact/(r13(ip)*rho(ip))
     630              :             sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
     631              :             sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
     632              :             sxxx = -280.0_dp/27.0_dp*s(ip)/(rho(ip)*rho(ip)*rho(ip))
     633              :             sxxy = 28.0_dp/9.0_dp*sfac*tact/(r13(ip)*rho(ip)*rho(ip)*rho(ip))
     634              :             e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + a3*fs(ip, 1) + 3.0_dp*a2*fs(ip, 2)*sx + &
     635              :                                 3.0_dp*a1*fs(ip, 3)*sx*sx + 3.0_dp*a1*fs(ip, 2)*sxx + &
     636              :                                 a0*fs(ip, 4)*sx*sx*sx + 3.0_dp*a0*fs(ip, 3)*sx*sxx + &
     637              :                                 a0*fs(ip, 2)*sxxx
     638              :             e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) + a2*fs(ip, 2)*sy + 2.0_dp*a1*fs(ip, 3)*sx*sy + &
     639              :                                   2.0_dp*a1*fs(ip, 2)*sxy + a0*fs(ip, 4)*sx*sx*sy + &
     640              :                                   2.0_dp*a0*fs(ip, 3)*sx*sxy + a0*fs(ip, 3)*sxx*sy + &
     641              :                                   a0*fs(ip, 2)*sxxy
     642              :             e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) + a1*fs(ip, 3)*sy*sy + a0*fs(ip, 4)*sx*sy*sy + &
     643              :                                     2.0_dp*a0*fs(ip, 3)*sxy*sy
     644              :             e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) + a0*fs(ip, 4)*sy*sy*sy
     645              : 
     646              :          END IF
     647              : 
     648              :       END DO
     649              : 
     650              : !$OMP     END PARALLEL DO
     651              : 
     652            0 :    END SUBROUTINE kex_p_3
     653              : 
     654              : ! Enhancement Factors
     655              : ! **************************************************************************************************
     656              : !> \brief ...
     657              : !> \param s ...
     658              : !> \param fs ...
     659              : !> \param m ...
     660              : ! **************************************************************************************************
     661            0 :    SUBROUTINE efactor_ol1(s, fs, m)
     662              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: s
     663              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: fs
     664              :       INTEGER, INTENT(IN)                                :: m
     665              : 
     666              :       INTEGER                                            :: ip
     667              :       REAL(KIND=dp)                                      :: t1, t2
     668              : 
     669            0 :       t1 = b*b/(72.0_dp*cf)
     670            0 :       t2 = 0.001878_dp*b
     671              : 
     672              : !$OMP     PARALLEL DO DEFAULT(NONE) &
     673              : !$OMP                 SHARED(s, m, fs, t1, t2) &
     674            0 : !$OMP                 PRIVATE(ip)
     675              : 
     676              :       DO ip = 1, SIZE(s)
     677              :          SELECT CASE (m)
     678              :          CASE (0)
     679              :             fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
     680              :          CASE (1)
     681              :             fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
     682              :             fs(ip, 2) = 2.0_dp*t1*s(ip) + t2
     683              :          CASE (2:3)
     684              :             fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
     685              :             fs(ip, 2) = 2.0_dp*t1*s(ip) + t2
     686              :             fs(ip, 3) = 2.0_dp*t1
     687              :          CASE DEFAULT
     688              :             CPABORT("Illegal order.")
     689              :          END SELECT
     690              :       END DO
     691              : 
     692              : !$OMP     END PARALLEL DO
     693              : 
     694            0 :       IF (m == 3) fs(:, 4) = 0.0_dp
     695              : 
     696            0 :    END SUBROUTINE efactor_ol1
     697              : ! **************************************************************************************************
     698              : !> \brief ...
     699              : !> \param s ...
     700              : !> \param fs ...
     701              : !> \param m ...
     702              : ! **************************************************************************************************
     703            0 :    SUBROUTINE efactor_ol2(s, fs, m)
     704              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: s
     705              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: fs
     706              :       INTEGER, INTENT(IN)                                :: m
     707              : 
     708              :       INTEGER                                            :: ip
     709              :       REAL(KIND=dp)                                      :: t1, t2, t3, y
     710              : 
     711            0 :       t1 = b*b/(72.0_dp*cf)
     712            0 :       t2 = 0.0245_dp*b
     713            0 :       t3 = 2.0_dp**f53*b
     714              : 
     715              : !$OMP     PARALLEL DO DEFAULT(NONE) &
     716              : !$OMP                 SHARED(s, m, t1, t2, t3, fs) &
     717            0 : !$OMP                 PRIVATE(ip,y)
     718              : 
     719              :       DO ip = 1, SIZE(s)
     720              :          y = 1.0_dp/(1.0_dp + t3*s(ip))
     721              :          SELECT CASE (m)
     722              :          CASE (0)
     723              :             fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
     724              :          CASE (1)
     725              :             fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
     726              :             fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
     727              :          CASE (2)
     728              :             fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
     729              :             fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
     730              :             fs(ip, 3) = 2.0_dp*(t1 - t2*t3*y*y*y)
     731              :          CASE (3)
     732              :             fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
     733              :             fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
     734              :             fs(ip, 3) = 2.0_dp*(t1 - t2*t3*y*y*y)
     735              :             fs(ip, 4) = 6.0_dp*t2*t3*t3*y*y*y*y
     736              :          CASE DEFAULT
     737              :             CPABORT("Illegal order.")
     738              :          END SELECT
     739              :       END DO
     740              : 
     741              : !$OMP     END PARALLEL DO
     742              : 
     743            0 :    END SUBROUTINE efactor_ol2
     744              : ! **************************************************************************************************
     745              : !> \brief ...
     746              : !> \param s ...
     747              : !> \param fs ...
     748              : !> \param m ...
     749              : ! **************************************************************************************************
     750          736 :    SUBROUTINE efactor_llp(s, fs, m)
     751              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: s
     752              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: fs
     753              :       INTEGER, INTENT(IN)                                :: m
     754              : 
     755              :       INTEGER                                            :: ip
     756              :       REAL(KIND=dp) :: as, bs, p, q, sas, sbs, t1, t10, t11, t12, t133, t16, t17, t19, t2, t20, &
     757              :          t22, t23, t26, t28, t29, t3, t30, t33, t36, t39, t4, t40, t42, t43, t45, t46, t47, t49, &
     758              :          t5, t50, t54, t55, t7, t71, t8, t9, x, ys
     759              : 
     760          736 :       p = 0.0044188_dp*b*b
     761          736 :       q = 0.0253_dp*b
     762              : 
     763              : !$OMP     PARALLEL DO DEFAULT(NONE) &
     764              : !$OMP                 SHARED(s, m, fs, p, q, b) &
     765              : !$OMP                 PRIVATE(ip,x,bs,sbs,as,sas,ys,t2,t4,t5,t8,t9,t10,t12) &
     766              : !$OMP                 PRIVATE(t1,t3,t7,t11,t16,t17,t20,t22,t23,t26,t30,t33) &
     767              : !$OMP                 PRIVATE(t40,t42,t43,t45,t46,t47,t49,t50,t71,t133) &
     768          736 : !$OMP                 PRIVATE(t19,t28,t29,t36,t39,t54,t55)
     769              : 
     770              :       DO ip = 1, SIZE(s)
     771              :          x = s(ip)
     772              :          bs = b*x
     773              :          sbs = SQRT(bs*bs + 1.0_dp)
     774              :          as = LOG(bs + sbs)
     775              :          sas = x*as
     776              :          ys = 1.0_dp/(1.0_dp + q*sas)
     777              :          SELECT CASE (m)
     778              :          CASE (0)
     779              :             fs(ip, 1) = 1.0_dp + p*x*x*ys
     780              :          CASE (1)
     781              :             fs(ip, 1) = 1.0_dp + p*x*x*ys
     782              :             t2 = q*x
     783              :             t4 = b**2
     784              :             t5 = x**2
     785              :             t8 = SQRT(1.0_dp + t4*t5)
     786              :             t9 = b*x + t8
     787              :             t10 = LOG(t9)
     788              :             t12 = 1.0_dp + t2*t10
     789              :             t17 = t12**2
     790              :             fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
     791              : 
     792              :          CASE (2)
     793              :             fs(ip, 1) = 1.0_dp + p*x*x*ys
     794              :             ! first der
     795              :             t2 = q*x
     796              :             t4 = b**2
     797              :             t5 = x**2
     798              :             t8 = SQRT(1.0_dp + t4*t5)
     799              :             t9 = b*x + t8
     800              :             t10 = LOG(t9)
     801              :             t12 = 1.0_dp + t2*t10
     802              :             t17 = t12**2
     803              :             fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
     804              : 
     805              :             ! second der
     806              :             t1 = q*x
     807              :             t3 = b**2
     808              :             t4 = x**2
     809              :             t7 = SQRT(1.0_dp + t3*t4)
     810              :             t8 = b*x + t7
     811              :             t9 = LOG(t8)
     812              :             t11 = 1.0_dp + t1*t9
     813              :             t16 = t11**2
     814              :             t17 = 1.0_dp/t16
     815              :             t20 = 1.0_dp/t7*t3
     816              :             t22 = b + t20*x
     817              :             t23 = 1/t8
     818              :             t26 = q*t9 + t1*t22*t23
     819              :             t30 = p*t4
     820              :             t33 = t26**2
     821              :             t40 = t7**2
     822              :             t43 = t3**2
     823              :             t49 = t22**2
     824              :             t50 = t8**2
     825              :             fs(ip, 3) = 2.0_dp*p/t11 - 4.0_dp*p*x*t17*t26 + 2.0_dp*t30/t16/ &
     826              :                         t11*t33 - t30*t17*(2.0_dp*q*t22*t23 + t1* &
     827              :                                            (-1.0_dp/t40/t7*t43*t4 + t20)*t23 - t1*t49/t50)
     828              : 
     829              :          CASE (3)
     830              : 
     831              :             fs(ip, 1) = 1.0_dp + p*x*x*ys
     832              :             ! first der
     833              :             t2 = q*x
     834              :             t4 = b**2
     835              :             t5 = x**2
     836              :             t8 = SQRT(1.0_dp + t4*t5)
     837              :             t9 = b*x + t8
     838              :             t10 = LOG(t9)
     839              :             t12 = 1.0_dp + t2*t10
     840              :             t17 = t12**2
     841              :             fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
     842              : 
     843              :             ! second der
     844              :             t1 = q*x
     845              :             t3 = b**2
     846              :             t4 = x**2
     847              :             t7 = SQRT(1.0_dp + t3*t4)
     848              :             t8 = b*x + t7
     849              :             t9 = LOG(t8)
     850              :             t11 = 1.0_dp + t1*t9
     851              :             t16 = t11**2
     852              :             t17 = 1.0_dp/t16
     853              :             t20 = 1.0_dp/t7*t3
     854              :             t22 = b + t20*x
     855              :             t23 = 1/t8
     856              :             t26 = q*t9 + t1*t22*t23
     857              :             t30 = p*t4
     858              :             t33 = t26**2
     859              :             t40 = t7**2
     860              :             t43 = t3**2
     861              :             t49 = t22**2
     862              :             t50 = t8**2
     863              :             fs(ip, 3) = 2.0_dp*p/t11 - 4.0_dp*p*x*t17*t26 + 2.0_dp*t30/t16/ &
     864              :                         t11*t33 - t30*t17*(2.0_dp*q*t22*t23 + t1* &
     865              :                                            (-1.0_dp/t40/t7*t43*t4 + t20)*t23 - t1*t49/t50)
     866              : 
     867              :             t1 = q*x
     868              :             t3 = b**2
     869              :             t4 = x**2
     870              :             t7 = SQRT(1 + t3*t4)
     871              :             t8 = b*x + t7
     872              :             t9 = LOG(t8)
     873              :             t11 = 1.0_dp + t1*t9
     874              :             t12 = t11**2
     875              :             t133 = 1.0_dp/t12
     876              :             t17 = 1.0_dp/t7*t3
     877              :             t19 = b + t17*x
     878              :             t20 = 1.0_dp/t8
     879              :             t23 = q*t9 + t1*t19*t20
     880              :             t26 = p*x
     881              :             t28 = 1.0_dp/t12/t11
     882              :             t29 = t23**2
     883              :             t36 = t7**2
     884              :             t39 = t3**2
     885              :             t40 = 1.0_dp/t36/t7*t39
     886              :             t42 = -t40*t4 + t17
     887              :             t45 = t19**2
     888              :             t46 = t8**2
     889              :             t47 = 1.0_dp/t46
     890              :             t50 = 2.0_dp*q*t19*t20 + t1*t42*t20 - t1*t45*t47
     891              :             t54 = p*t4
     892              :             t55 = t12**2
     893              :             t71 = t36**2
     894              :             fs(ip, 4) = &
     895              :                -6.0_dp*p*t133*t23 + 12.0_dp*t26*t28*t29 - &
     896              :                6.0_dp*t26*t133*t50 - 6.0_dp*t54/t55*t29*t23 + &
     897              :                6.0_dp*t54*t28*t23*t50 - t54*t133* &
     898              :                (3.0_dp*q*t42*t20 - 3.0_dp*q*t45*t47 + 3.0_dp*t1* &
     899              :                 (1.0_dp/t71/t7*t39*t3*t4*x - t40*x)*t20 - &
     900              :                 3.0_dp*t1*t42*t47*t19 + 2.0_dp*t1*t45*t19/t46/t8)
     901              : 
     902              :          CASE DEFAULT
     903              :             CPABORT("Illegal order.")
     904              :          END SELECT
     905              :       END DO
     906              : 
     907              : !$OMP     END PARALLEL DO
     908              : 
     909          736 :    END SUBROUTINE efactor_llp
     910              : ! **************************************************************************************************
     911              : !> \brief ...
     912              : !> \param s ...
     913              : !> \param fs ...
     914              : !> \param m ...
     915              : !> \param f2_lsd ...
     916              : ! **************************************************************************************************
     917           54 :    SUBROUTINE efactor_pw86(s, fs, m, f2_lsd)
     918              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: s
     919              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: fs
     920              :       INTEGER, INTENT(IN)                                :: m
     921              :       REAL(dp), OPTIONAL                                 :: f2_lsd
     922              : 
     923              :       INTEGER                                            :: ip
     924              :       REAL(KIND=dp)                                      :: f15, ff, p0, p1, p15, p2, p3, s1, s2, &
     925              :                                                             s4, s6, t1, t2, t3
     926              : 
     927           54 :       t1 = 1.296_dp
     928           54 :       t2 = 14.0_dp
     929           54 :       t3 = 0.2_dp
     930           54 :       f15 = 1.0_dp/15.0_dp
     931           54 :       ff = 1.0_dp
     932           54 :       IF (PRESENT(f2_lsd)) ff = f2_lsd
     933              : 
     934              : !$OMP     PARALLEL DO DEFAULT(NONE) &
     935              : !$OMP                 SHARED(s, fs, m, t1, t2, t3, f15, ff) &
     936           54 : !$OMP                 PRIVATE(ip, s1, s2, s4, s6, p0, p1, p2, p3, p15)
     937              : 
     938              :       DO ip = 1, SIZE(s)
     939              :          s1 = s(ip)*ff
     940              :          s2 = s1*s1
     941              :          s4 = s2*s2
     942              :          s6 = s2*s4
     943              :          SELECT CASE (m)
     944              :          CASE (0)
     945              :             p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
     946              :             fs(ip, 1) = p0**f15
     947              :          CASE (1)
     948              :             p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
     949              :             p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
     950              :             p15 = p0**f15
     951              :             fs(ip, 1) = p15
     952              :             fs(ip, 2) = f15*p1*p15/p0
     953              :          CASE (2)
     954              :             p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
     955              :             p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
     956              :             p2 = ff*ff*(2.0_dp*t1 + 12.0_dp*t2*s2 + 30.0_dp*t3*s4)
     957              :             p15 = p0**f15
     958              :             fs(ip, 1) = p15
     959              :             fs(ip, 2) = f15*p1*p15/p0
     960              :             fs(ip, 3) = f15*p15/p0*(p2 - 14.0_dp/15.0_dp*p1*p1/p0)
     961              :          CASE (3)
     962              :             p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
     963              :             p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
     964              :             p2 = ff*ff*(2.0_dp*t1 + 12.0_dp*t2*s2 + 30.0_dp*t3*s4)
     965              :             p3 = s1*ff*ff*ff*(24.0_dp*t2 + 120.0_dp*t3*s2)
     966              :             p15 = p0**f15
     967              :             fs(ip, 1) = p15
     968              :             fs(ip, 2) = f15*p1*p15/p0
     969              :             fs(ip, 3) = f15*p15/p0*(p2 - 14.0_dp/15.0_dp*p1*p1/p0)
     970              :             fs(ip, 4) = f15*p15/p0*(-14.0_dp*f15*p1*p1/p0 + 14.0_dp*14.0_dp*f15*p1*p1*p1/p0/p0 + &
     971              :                                     p3 - 14.0_dp*p2*p1/p0 + 14.0_dp*p1*p1/p0/p0)
     972              :          CASE DEFAULT
     973              :             CPABORT("Illegal order.")
     974              :          END SELECT
     975              :       END DO
     976              : 
     977              : !$OMP     END PARALLEL DO
     978              : 
     979           54 :    END SUBROUTINE efactor_pw86
     980              : ! **************************************************************************************************
     981              : !> \brief ...
     982              : !> \param s ...
     983              : !> \param fs ...
     984              : !> \param m ...
     985              : ! **************************************************************************************************
     986          174 :    SUBROUTINE efactor_t92(s, fs, m)
     987              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: s
     988              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: fs
     989              :       INTEGER, INTENT(IN)                                :: m
     990              : 
     991              :       INTEGER                                            :: ip
     992              :       REAL(KIND=dp)                                      :: a1, a2, as, asp, asp2, asp3, bs, p, q, &
     993              :                                                             s1, s2, sas, sbs, sbs3, sbs5, t0, w1, &
     994              :                                                             x, ys
     995              : 
     996          174 :       p = 0.0055_dp*b*b
     997          174 :       q = 0.0253_dp*b
     998          174 :       a1 = 0.072_dp*b
     999          174 :       a2 = 2.0_dp**f53*b
    1000              : 
    1001              : !$OMP     PARALLEL DO DEFAULT(NONE) &
    1002              : !$OMP                 SHARED(s, fs, m, p, q, a1, a2, b) &
    1003              : !$OMP                 PRIVATE(ip, x, bs, sbs, sas, ys, asp, sbs3, asp2, sbs5) &
    1004          174 : !$OMP                 PRIVATE(asp3, as, s2, s1, t0, w1)
    1005              : 
    1006              :       DO ip = 1, SIZE(s)
    1007              :          x = s(ip)
    1008              :          bs = b*x
    1009              :          sbs = SQRT(bs*bs + 1.0_dp)
    1010              :          as = LOG(bs + sbs)
    1011              :          sas = x*as
    1012              :          ys = 1.0_dp/(1.0_dp + q*sas)
    1013              :          SELECT CASE (m)
    1014              :          CASE (0)
    1015              :             fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
    1016              :          CASE (1)
    1017              :             asp = as + bs/sbs
    1018              :             fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
    1019              :             fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
    1020              :          CASE (2)
    1021              :             asp = as + bs/sbs
    1022              :             sbs3 = sbs*sbs*sbs
    1023              :             asp2 = 2.0_dp*b/sbs - b*bs*bs/sbs3
    1024              :             fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
    1025              :             fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
    1026              :             fs(ip, 3) = 2.0_dp*p*ys - p*q*x*(4.0_dp*asp + x*asp2)*ys*ys + &
    1027              :                         2.0_dp*p*q*q*x*x*asp*asp*ys*ys*ys + 2.0_dp*a1*a2/(1 + a2*x)**3
    1028              :          CASE (3)
    1029              :             asp = as + bs/sbs
    1030              :             sbs3 = sbs*sbs*sbs
    1031              :             sbs5 = sbs3*sbs*sbs
    1032              :             asp2 = 2.0_dp*b/sbs - b*bs*bs/sbs3
    1033              :             asp3 = -4.0_dp*b*b*bs/sbs3 + 3.0_dp*b*b*bs*bs*bs/sbs5
    1034              :             w1 = (4.0_dp*asp + x*asp2)
    1035              :             fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
    1036              :             fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
    1037              :             fs(ip, 3) = 2.0_dp*p*ys - p*q*x*w1*ys*ys + &
    1038              :                         2.0_dp*p*q*q*x*x*asp*asp*ys*ys*ys + 2.0_dp*a1*a2/(1 + a2*x)**3
    1039              : 
    1040              :             s2 = -6*p/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**2*(q*LOG(b*x + SQRT(1 + b**2*x**2)) + &
    1041              :                                                      q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/(b*x + SQRT(1 + b**2*x**2))) + 12*p*x/ &
    1042              :                  (1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**3*(q*LOG(b*x + SQRT(1 + b**2*x**2)) + &
    1043              :                                                               q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/(b*x + SQRT(1 + b**2*x**2)))**2
    1044              :             s1 = s2 - 6*p*x/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**2*(2*q*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/ &
    1045              :                              (b*x + SQRT(1 + b**2*x**2)) + q*x*(-1/SQRT(1 + b**2*x**2)**3*b**4*x**2 + 1/SQRT(1 + b**2*x**2)*b**2)/ &
    1046              :                                                           (b*x + SQRT(1 + b**2*x**2)) - q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)**2/ &
    1047              :                                             (b*x + SQRT(1 + b**2*x**2))**2) - 6*p*x**2/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**4 &
    1048              :                  *(q*LOG(b*x + SQRT(1 + b**2*x**2)) + q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/(b*x + SQRT(1 + b**2*x**2)))**3
    1049              :             s2 = s1 + 6*p*x**2/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**3*(q*LOG(b*x + SQRT(1 + b**2*x**2)) + &
    1050              :                        q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/(b*x + SQRT(1 + b**2*x**2)))*(2*q*(b + 1/SQRT(1 + b**2*x**2)*b**2*x) &
    1051              :                                   /(b*x + SQRT(1 + b**2*x**2)) + q*x*(-1/SQRT(1 + b**2*x**2)**3*b**4*x**2 + 1/SQRT(1 + b**2*x**2)* &
    1052              :                        b**2)/(b*x + SQRT(1 + b**2*x**2)) - q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)**2/(b*x + SQRT(1 + b**2*x**2))**2)
    1053              :             t0 = s2 - p*x**2/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**2*(3*q*(-1/SQRT(1 + b**2*x**2)**3*b**4*x**2 + &
    1054              :                               1/SQRT(1 + b**2*x**2)*b**2)/(b*x + SQRT(1 + b**2*x**2)) - 3*q*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)**2/ &
    1055              :                       (b*x + SQRT(1 + b**2*x**2))**2 + q*x*(3/SQRT(1 + b**2*x**2)**5*b**6*x**3 - 3/SQRT(1 + b**2*x**2)**3*b**4*x)/ &
    1056              :                            (b*x + SQRT(1 + b**2*x**2)) - 3*q*x*(-1/SQRT(1 + b**2*x**2)**3*b**4*x**2 + 1/SQRT(1 + b**2*x**2)*b**2)/ &
    1057              :                              (b*x + SQRT(1 + b**2*x**2))**2*(b + 1/SQRT(1 + b**2*x**2)*b**2*x) + 2*q*x*(b + 1/SQRT(1 + b**2*x**2)* &
    1058              :                                   b**2*x)**3/(b*x + SQRT(1 + b**2*x**2))**3) - 6*a1/(1 + a2*x)**3*a2**2 + 6*a1*x/(1 + a2*x)**4*a2**3
    1059              : 
    1060              :             fs(ip, 4) = t0
    1061              : 
    1062              :          CASE DEFAULT
    1063              :             CPABORT("Illegal order")
    1064              :          END SELECT
    1065              :       END DO
    1066              : 
    1067              : !$OMP     END PARALLEL DO
    1068              : 
    1069          174 :    END SUBROUTINE efactor_t92
    1070              : ! **************************************************************************************************
    1071              : !> \brief ...
    1072              : !> \param s ...
    1073              : !> \param fs ...
    1074              : !> \param m ...
    1075              : !> \param pset ...
    1076              : !> \param f2_lsd ...
    1077              : ! **************************************************************************************************
    1078          576 :    SUBROUTINE efactor_pbex(s, fs, m, pset, f2_lsd)
    1079              : 
    1080              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: s
    1081              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: fs
    1082              :       INTEGER, INTENT(IN)                                :: m, pset
    1083              :       REAL(dp), OPTIONAL                                 :: f2_lsd
    1084              : 
    1085              :       REAL(KIND=dp), PARAMETER                           :: kappa1 = 0.804_dp, kappa2 = 1.245_dp, &
    1086              :                                                             mu = 0.2195149727645171_dp
    1087              : 
    1088              :       INTEGER                                            :: ip
    1089              :       REAL(KIND=dp)                                      :: f0, mk, x, x2, y
    1090              : 
    1091              :       IF (pset == 1) mk = mu/kappa1
    1092          576 :       IF (pset == 2) mk = mu/kappa2
    1093              : 
    1094          576 :       f0 = 1.0_dp/tact
    1095          576 :       IF (PRESENT(f2_lsd)) f0 = f2_lsd
    1096              : 
    1097              : !$OMP     PARALLEL DO DEFAULT(NONE) &
    1098              : !$OMP                 SHARED(s, m, fs, f0, mk) &
    1099          576 : !$OMP                 PRIVATE(ip,x,x2,y)
    1100              : 
    1101              :       DO ip = 1, SIZE(s)
    1102              :          x = s(ip)*f0
    1103              :          x2 = x*x
    1104              :          y = 1.0_dp/(1.0_dp + mk*x2)
    1105              :          SELECT CASE (m)
    1106              :          CASE (0)
    1107              :             fs(ip, 1) = 1.0_dp + mu*x2*y
    1108              :          CASE (1)
    1109              :             fs(ip, 1) = 1.0_dp + mu*x2*y
    1110              :             fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
    1111              :          CASE (2)
    1112              :             fs(ip, 1) = 1.0_dp + mu*x2*y
    1113              :             fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
    1114              :             fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
    1115              :          CASE (3)
    1116              :             fs(ip, 1) = 1.0_dp + mu*x2*y
    1117              :             fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
    1118              :             fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
    1119              :             fs(ip, 4) = 24.0_dp*mu*mk*x*(mk*x2 - 1.0_dp)*y*y*y*y*f0*f0*f0
    1120              :          CASE DEFAULT
    1121              :             CPABORT("Illegal order.")
    1122              :          END SELECT
    1123              :       END DO
    1124              : 
    1125              : !$OMP     END PARALLEL DO
    1126              : 
    1127          576 :    END SUBROUTINE efactor_pbex
    1128              : 
    1129              : ! **************************************************************************************************
    1130              : !> \brief ...
    1131              : !> \param s ...
    1132              : !> \param fs ...
    1133              : !> \param m ...
    1134              : !> \param pset ...
    1135              : !> \param f2_lsd ...
    1136              : ! **************************************************************************************************
    1137          108 :    SUBROUTINE efactor_pw91(s, fs, m, pset, f2_lsd)
    1138              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: s
    1139              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: fs
    1140              :       INTEGER, INTENT(IN)                                :: m, pset
    1141              :       REAL(dp), OPTIONAL                                 :: f2_lsd
    1142              : 
    1143              :       INTEGER                                            :: ip
    1144              :       REAL(dp) :: ff, t1, t10, t101, t106, t109, t111, t113, t119, t12, t123, t124, t13, t14, t15, &
    1145              :          t16, t17, t18, t19, t2, t20, t21, t22, t23, t25, t26, t27, t28, t29, t3, t30, t31, t33, &
    1146              :          t35, t37, t38, t39, t4, t40, t44, t47, t48, t5, t50, t51, t53, t55, t56, t57, t58, t59, &
    1147              :          t6, t60, t64, t65, t69, t7, t70, t71, t73, t77, t78, t8, t80, t82, t9, t90, t93, t94, &
    1148              :          t96, t98
    1149              :       REAL(KIND=dp)                                      :: a1, a2, a3, a4, a5, bb, o, pa(6, 2), x
    1150              : 
    1151              : ! parameter set 1: Perdew-Wang
    1152              : ! parameter set 2: Lembarki-Chermette
    1153              : 
    1154              :       pa(1:6, 1) = (/0.19645_dp, 0.2743_dp, &
    1155              :                      0.1508_dp, 100.0_dp, &
    1156          756 :                      7.7956_dp, 0.004_dp/)
    1157              :       pa(1:6, 2) = (/0.093907_dp, 0.26608_dp, &
    1158              :                      0.0809615_dp, 100.0_dp, &
    1159          756 :                      76.320_dp, 0.57767e-4_dp/)
    1160          108 :       o = 1.0_dp
    1161          108 :       ff = 1.0_dp
    1162          108 :       IF (PRESENT(f2_lsd)) ff = f2_lsd
    1163          108 :       a1 = pa(1, pset)*FF
    1164          108 :       a2 = pa(2, pset)*ff*ff
    1165          108 :       a3 = pa(3, pset)*ff*ff
    1166          108 :       a4 = pa(4, pset)*ff*ff
    1167          108 :       bb = pa(5, pset)*ff
    1168              : !   it should be valid also for lsd
    1169          108 :       a5 = pa(6, pset)*ff*ff*ff*ff
    1170              : 
    1171              : !$OMP     PARALLEL DEFAULT(NONE) &
    1172              : !$OMP              SHARED(s, fs, m, a1, a2, a3, a4, a5, bb, pa, o, ff) &
    1173              : !$OMP              PRIVATE(x,t1,t10,t101,t106,t109,t111,t113,t119,t12,t123) &
    1174              : !$OMP              PRIVATE(t124,t13,t14,t15,t16,t17,t18,t19,t2,t20,t21,t22) &
    1175              : !$OMP              PRIVATE(t23,t25,t26,t27,t28,t29,t3,t30,t31,t33,t35,t37) &
    1176              : !$OMP              PRIVATE(t38,t39,t4,t40,t44,t47,t48,t5,t50,t51,t53,t55) &
    1177              : !$OMP              PRIVATE(t56,t57,t58,t59,t6,t60,t64,t65,t69,t7,t70,t71) &
    1178          108 : !$OMP              PRIVATE(t73,t77,t78,t8,t80,t82,t9,t90,t93,t94,t96,t98,ip)
    1179              : 
    1180              :       IF (m >= 0) THEN
    1181              : !$OMP       DO
    1182              :          DO ip = 1, SIZE(s)
    1183              :             x = s(ip)
    1184              :             t3 = bb**2
    1185              :             t4 = x**2
    1186              :             t7 = SQRT(o + t3*t4)
    1187              :             t9 = LOG(bb*x + t7)
    1188              :             t10 = a1*x*t9
    1189              :             t12 = EXP(-a4*t4)
    1190              :             t17 = t4**2
    1191              :             fs(ip, 1) = (o + t10 + (a2 - a3*t12)*t4)/(o + t10 + a5*t17)
    1192              :          END DO
    1193              : !$OMP       END DO
    1194              :       END IF
    1195              :       IF (m >= 1) THEN
    1196              : !$OMP       DO
    1197              :          DO ip = 1, SIZE(s)
    1198              :             x = s(ip)
    1199              :             t2 = bb**2
    1200              :             t3 = x**2
    1201              :             t6 = SQRT(o + t2*t3)
    1202              :             t7 = bb*x + t6
    1203              :             t8 = LOG(t7)
    1204              :             t9 = a1*t8
    1205              :             t10 = a1*x
    1206              :             t17 = t10*(bb + 1/t6*t2*x)/t7
    1207              :             t19 = t3*x
    1208              :             t21 = EXP(-a4*t3)
    1209              :             t26 = a2 - a3*t21
    1210              :             t30 = t10*t8
    1211              :             t31 = t3**2
    1212              :             t33 = o + t30 + a5*t31
    1213              :             t38 = t33**2
    1214              :             fs(ip, 2) = &
    1215              :                (t9 + t17 + 2._dp*a3*a4*t19*t21 + 2._dp*t26*x)/ &
    1216              :                t33 - (o + t30 + t26*t3)/t38*(t9 + t17 + 4._dp*a5*t19)
    1217              :          END DO
    1218              : !$OMP       END DO
    1219              :       END IF
    1220              :       IF (m >= 2) THEN
    1221              : !$OMP       DO
    1222              :          DO ip = 1, SIZE(s)
    1223              :             x = s(ip)
    1224              :             t1 = bb**2
    1225              :             t2 = x**2
    1226              :             t5 = SQRT(o + t1*t2)
    1227              :             t7 = o/t5*t1
    1228              :             t9 = bb + t7*x
    1229              :             t12 = bb*x + t5
    1230              :             t13 = o/t12
    1231              :             t15 = 2._dp*a1*t9*t13
    1232              :             t16 = a1*x
    1233              :             t17 = t5**2
    1234              :             t20 = t1**2
    1235              :             t25 = t16*(-o/t17/t5*t20*t2 + t7)*t13
    1236              :             t26 = t9**2
    1237              :             t27 = t12**2
    1238              :             t30 = t16*t26/t27
    1239              :             t31 = a3*a4
    1240              :             t33 = EXP(-a4*t2)
    1241              :             t37 = a4**2
    1242              :             t39 = t2**2
    1243              :             t44 = a3*t33
    1244              :             t47 = LOG(t12)
    1245              :             t48 = t16*t47
    1246              :             t50 = o + t48 + a5*t39
    1247              :             t53 = a1*t47
    1248              :             t55 = t16*t9*t13
    1249              :             t56 = t2*x
    1250              :             t60 = a2 - t44
    1251              :             t64 = t50**2
    1252              :             t65 = o/t64
    1253              :             t69 = t53 + t55 + 4._dp*a5*t56
    1254              :             t73 = o + t48 + t60*t2
    1255              :             t77 = t69**2
    1256              :             fs(ip, 3) = &
    1257              :                (t15 + t25 - t30 + 10._dp*t31*t2*t33 - 4._dp*a3*t37*t39*t33 + &
    1258              :                 2._dp*a2 - 2._dp*t44)/t50 - 2._dp* &
    1259              :                (t53 + t55 + 2._dp*t31*t56*t33 + 2._dp*t60*x)* &
    1260              :                t65*t69 + 2._dp*t73/t64/t50*t77 - t73*t65*(t15 + t25 - t30 + 12._dp*a5*t2)
    1261              :          END DO
    1262              : !$OMP       END DO
    1263              :       END IF
    1264              :       IF (m >= 3) THEN
    1265              : !$OMP       DO
    1266              :          DO ip = 1, SIZE(s)
    1267              :             x = s(ip)
    1268              :             t1 = bb**2
    1269              :             t2 = x**2
    1270              :             t5 = SQRT(0.1e1_dp + t1*t2)
    1271              :             t6 = t5**2
    1272              :             t9 = t1**2
    1273              :             t10 = 1/t6/t5*t9
    1274              :             t13 = 1/t5*t1
    1275              :             t14 = -t10*t2 + t13
    1276              :             t17 = bb*x + t5
    1277              :             t18 = 1/t17
    1278              :             t20 = 3*a1*t14*t18
    1279              :             t22 = bb + t13*x
    1280              :             t23 = t22**2
    1281              :             t25 = t17**2
    1282              :             t26 = 1/t25
    1283              :             t28 = 3*a1*t23*t26
    1284              :             t29 = a1*x
    1285              :             t30 = t6**2
    1286              :             t35 = t2*x
    1287              :             t40 = 3*t29*(1/t30/t5*t1*t9*t35 - t10*x)*t18
    1288              :             t44 = 3*t29*t14*t26*t22
    1289              :             t50 = 2*t29*t23*t22/t25/t17
    1290              :             t51 = a3*a4
    1291              :             t53 = EXP(-a4*t2)
    1292              :             t57 = a4**2
    1293              :             t58 = a3*t57
    1294              :             t59 = t35*t53
    1295              :             t64 = t2**2
    1296              :             t70 = LOG(t17)
    1297              :             t71 = t29*t70
    1298              :             t73 = 0.1e1_dp + t71 + a5*t64
    1299              :             t78 = 2*a1*t22*t18
    1300              :             t80 = t29*t14*t18
    1301              :             t82 = t29*t23*t26
    1302              :             t90 = a3*t53
    1303              :             t93 = t73**2
    1304              :             t94 = 1/t93
    1305              :             t96 = a1*t70
    1306              :             t98 = t29*t18*t22
    1307              :             t101 = t96 + t98 + 4*a5*t35
    1308              :             t106 = a2 - t90
    1309              :             t109 = t96 + t98 + 2*t51*t59 + 2*t106*x
    1310              :             t111 = 1/t93/t73
    1311              :             t113 = t101**2
    1312              :             t119 = t78 + t80 - t82 + 12*a5*t2
    1313              :             t123 = 0.1e1_dp + t71 + t106*t2
    1314              :             t124 = t93**2
    1315              :             fs(ip, 4) = &
    1316              :                (t20 - t28 + t40 - t44 + t50 + 24*t51*x*t53 - 36._dp*t58*t59 + 8._dp*a3*t57*a4*t64* &
    1317              :                 x*t53)/t73 - 3._dp*(t78 + t80 - t82 + 10._dp*t51*t2*t53 - &
    1318              :                                     4._dp*t58*t64*t53 + 2._dp*a2 - 2._dp*t90)*t94*t101 + &
    1319              :                6._dp*t109*t111*t113 - 3._dp*t109*t94*t119 - 6*t123/t124*t113*t101 + &
    1320              :                6._dp*t123*t111*t101*t119 - t123*t94*(t20 - t28 + t40 - t44 + t50 + 24._dp*a5*x)
    1321              :          END DO
    1322              : !$OMP       END DO
    1323              :       END IF
    1324              : 
    1325              : !$OMP     END PARALLEL
    1326              : 
    1327          108 :    END SUBROUTINE efactor_pw91
    1328              : 
    1329              : END MODULE xc_ke_gga
    1330              : 
        

Generated by: LCOV version 2.0-1