LCOV - code coverage report
Current view: top level - src/xc - xc_pade.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 63.5 % 244 155
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 LDA functional in the Pade approximation
      10              : !>      Literature: S. Goedecker, M. Teter and J. Hutter,
      11              : !>                  Phys. Rev. B 54, 1703 (1996)
      12              : !> \note
      13              : !>      Order of derivatives is: LDA 0; 1; 2; 3;
      14              : !>                               LSD 0; a  b; aa ab bb; aaa aab abb bbb;
      15              : !> \par History
      16              : !>      JGH (26.02.2003) : OpenMP enabled
      17              : !> \author JGH (15.02.2002)
      18              : ! **************************************************************************************************
      19              : MODULE xc_pade
      20              :    USE bibliography,                    ONLY: Goedecker1996,&
      21              :                                               cite_reference
      22              :    USE kinds,                           ONLY: dp
      23              :    USE pw_types,                        ONLY: pw_r3d_rs_type
      24              :    USE xc_derivative_desc,              ONLY: 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_fx,&
      32              :                                               calc_rs,&
      33              :                                               calc_rs_pw,&
      34              :                                               set_util
      35              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      36              :    USE xc_rho_set_types,                ONLY: xc_rho_set_type
      37              : #include "../base/base_uses.f90"
      38              : 
      39              :    IMPLICIT NONE
      40              : 
      41              :    PRIVATE
      42              : 
      43              :    REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
      44              :                                f23 = 2.0_dp*f13, &
      45              :                                f43 = 4.0_dp*f13
      46              : 
      47              :    REAL(KIND=dp), PARAMETER :: a0 = 0.4581652932831429E+0_dp, &
      48              :                                a1 = 0.2217058676663745E+1_dp, &
      49              :                                a2 = 0.7405551735357053E+0_dp, &
      50              :                                a3 = 0.1968227878617998E-1_dp, &
      51              :                                b1 = 1.0000000000000000E+0_dp, &
      52              :                                b2 = 0.4504130959426697E+1_dp, &
      53              :                                b3 = 0.1110667363742916E+1_dp, &
      54              :                                b4 = 0.2359291751427506E-1_dp
      55              : 
      56              :    REAL(KIND=dp), PARAMETER :: da0 = 0.119086804055547E+0_dp, &
      57              :                                da1 = 0.6157402568883345E+0_dp, &
      58              :                                da2 = 0.1574201515892867E+0_dp, &
      59              :                                da3 = 0.3532336663397157E-2_dp, &
      60              :                                db1 = 0.0000000000000000E+0_dp, &
      61              :                                db2 = 0.2673612973836267E+0_dp, &
      62              :                                db3 = 0.2052004607777787E+0_dp, &
      63              :                                db4 = 0.4200005045691381E-2_dp
      64              : 
      65              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pade'
      66              : 
      67              :    PUBLIC :: pade_lda_pw_eval, pade_lsd_pw_eval, pade_info, pade_init, pade_fxc_eval
      68              : 
      69              :    REAL(KIND=dp) :: eps_rho
      70              :    LOGICAL :: debug_flag
      71              : 
      72              : CONTAINS
      73              : 
      74              : ! **************************************************************************************************
      75              : !> \brief ...
      76              : !> \param cutoff ...
      77              : !> \param debug ...
      78              : ! **************************************************************************************************
      79        78656 :    SUBROUTINE pade_init(cutoff, debug)
      80              : 
      81              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
      82              :       LOGICAL, INTENT(IN), OPTIONAL                      :: debug
      83              : 
      84        78656 :       eps_rho = cutoff
      85        78656 :       CALL set_util(cutoff)
      86              : 
      87        78656 :       CALL cite_reference(Goedecker1996)
      88              : 
      89        78656 :       IF (PRESENT(debug)) THEN
      90            0 :          debug_flag = debug
      91              :       ELSE
      92        78656 :          debug_flag = .FALSE.
      93              :       END IF
      94              : 
      95        78656 :    END SUBROUTINE pade_init
      96              : 
      97              : ! **************************************************************************************************
      98              : !> \brief ...
      99              : !> \param reference ...
     100              : !> \param shortform ...
     101              : !> \param lsd ...
     102              : !> \param needs ...
     103              : !> \param max_deriv ...
     104              : ! **************************************************************************************************
     105        76759 :    SUBROUTINE pade_info(reference, shortform, lsd, needs, max_deriv)
     106              : 
     107              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
     108              :       LOGICAL, INTENT(IN), OPTIONAL                      :: lsd
     109              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
     110              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
     111              : 
     112        76759 :       IF (PRESENT(reference)) THEN
     113              :          reference = "S. Goedecker, M. Teter and J. Hutter," &
     114          517 :                      //" Phys. Rev. B 54, 1703 (1996)"
     115              :       END IF
     116        76759 :       IF (PRESENT(shortform)) THEN
     117          517 :          shortform = "S. Goedecker et al., PRB 54, 1703 (1996)"
     118              :       END IF
     119              : 
     120        76759 :       IF (PRESENT(needs)) THEN
     121        76242 :          IF (.NOT. PRESENT(lsd)) &
     122            0 :             CPABORT("Arguments mismatch.")
     123        76242 :          IF (lsd) THEN
     124        12739 :             needs%rho_spin = .TRUE.
     125              :          ELSE
     126        63503 :             needs%rho = .TRUE.
     127              :          END IF
     128              :       END IF
     129              : 
     130        76759 :       IF (PRESENT(max_deriv)) max_deriv = 3
     131              : 
     132        76759 :    END SUBROUTINE pade_info
     133              : 
     134              : ! **************************************************************************************************
     135              : !> \brief ...
     136              : !> \param deriv_set ...
     137              : !> \param rho_set ...
     138              : !> \param order ...
     139              : ! **************************************************************************************************
     140        65691 :    SUBROUTINE pade_lda_pw_eval(deriv_set, rho_set, order)
     141              : 
     142              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     143              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     144              :       INTEGER, INTENT(IN), OPTIONAL                      :: order
     145              : 
     146              :       INTEGER                                            :: n
     147              :       LOGICAL                                            :: calc(0:4)
     148        65691 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: rs
     149              :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     150        65691 :          POINTER                                         :: e_0, e_r, e_rr, e_rrr
     151              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     152              : 
     153        65691 :       calc = .FALSE.
     154       199715 :       IF (order >= 0) calc(0:order) = .TRUE.
     155        65691 :       IF (order < 0) calc(-order) = .TRUE.
     156              : 
     157       262764 :       n = PRODUCT(rho_set%local_bounds(2, :) - rho_set%local_bounds(1, :) + (/1, 1, 1/))
     158       197073 :       ALLOCATE (rs(n))
     159              : 
     160        65691 :       CALL calc_rs_pw(rho_set%rho, rs, n)
     161        65691 :       IF (calc(0) .AND. calc(1)) THEN
     162              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     163        62727 :                                          allocate_deriv=.TRUE.)
     164        62727 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     165              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     166        62727 :                                          allocate_deriv=.TRUE.)
     167        62727 :          CALL xc_derivative_get(deriv, deriv_data=e_r)
     168        62727 :          CALL pade_lda_01(n, rho_set%rho, rs, e_0, e_r)
     169         2964 :       ELSE IF (calc(0)) THEN
     170              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     171         2964 :                                          allocate_deriv=.TRUE.)
     172         2964 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     173         2964 :          CALL pade_lda_0(n, rho_set%rho, rs, e_0)
     174            0 :       ELSE IF (calc(1)) THEN
     175              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     176            0 :                                          allocate_deriv=.TRUE.)
     177            0 :          CALL xc_derivative_get(deriv, deriv_data=e_r)
     178            0 :          CALL pade_lda_1(n, rho_set%rho, rs, e_r)
     179              :       END IF
     180        65691 :       IF (calc(2)) THEN
     181              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     182         5606 :                                          allocate_deriv=.TRUE.)
     183         5606 :          CALL xc_derivative_get(deriv, deriv_data=e_rr)
     184         5606 :          CALL pade_lda_2(n, rho_set%rho, rs, e_rr)
     185              :       END IF
     186        65691 :       IF (calc(3)) THEN
     187              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
     188            0 :                                          allocate_deriv=.TRUE.)
     189            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rrr)
     190            0 :          CALL pade_lda_3(n, rho_set%rho, rs, e_rrr)
     191              :       END IF
     192              : 
     193        65691 :       DEALLOCATE (rs)
     194              : 
     195        65691 :    END SUBROUTINE pade_lda_pw_eval
     196              : 
     197              : ! **************************************************************************************************
     198              : !> \brief ...
     199              : !> \param deriv_set ...
     200              : !> \param rho_set ...
     201              : !> \param order ...
     202              : ! **************************************************************************************************
     203        12963 :    SUBROUTINE pade_lsd_pw_eval(deriv_set, rho_set, order)
     204              : 
     205              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     206              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     207              :       INTEGER, INTENT(IN), OPTIONAL                      :: order
     208              : 
     209              :       INTEGER                                            :: i, j, k
     210              :       LOGICAL                                            :: calc(0:4)
     211              :       REAL(KIND=dp)                                      :: rhoa, rhob, rs
     212              :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     213        12963 :          POINTER                                         :: e_0, e_ra, e_rara, e_rarara, e_rararb, &
     214        12963 :                                                             e_rarb, e_rarbrb, e_rb, e_rbrb, &
     215        12963 :                                                             e_rbrbrb
     216              :       REAL(KIND=dp), DIMENSION(4)                        :: fx
     217              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     218              : 
     219        12963 :       calc = .FALSE.
     220        38199 :       IF (order >= 0) calc(0:order) = .TRUE.
     221        12963 :       IF (order < 0) calc(-order) = .TRUE.
     222              : 
     223        12963 :       IF (calc(0)) THEN
     224              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     225        12963 :                                          allocate_deriv=.TRUE.)
     226        12963 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     227              :       END IF
     228        12963 :       IF (calc(1)) THEN
     229              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
     230        11537 :                                          allocate_deriv=.TRUE.)
     231        11537 :          CALL xc_derivative_get(deriv, deriv_data=e_ra)
     232              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
     233        11537 :                                          allocate_deriv=.TRUE.)
     234        11537 :          CALL xc_derivative_get(deriv, deriv_data=e_rb)
     235              :       END IF
     236        12963 :       IF (calc(2)) THEN
     237              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
     238          736 :                                          allocate_deriv=.TRUE.)
     239          736 :          CALL xc_derivative_get(deriv, deriv_data=e_rara)
     240              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
     241          736 :                                          allocate_deriv=.TRUE.)
     242          736 :          CALL xc_derivative_get(deriv, deriv_data=e_rarb)
     243              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
     244          736 :                                          allocate_deriv=.TRUE.)
     245          736 :          CALL xc_derivative_get(deriv, deriv_data=e_rbrb)
     246              :       END IF
     247        12963 :       IF (calc(3)) THEN
     248              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
     249            0 :                                          allocate_deriv=.TRUE.)
     250            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rarara)
     251              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
     252            0 :                                          allocate_deriv=.TRUE.)
     253            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rararb)
     254              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
     255            0 :                                          allocate_deriv=.TRUE.)
     256            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rarbrb)
     257              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
     258            0 :                                          allocate_deriv=.TRUE.)
     259            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rbrbrb)
     260              :       END IF
     261              : 
     262              : !$OMP PARALLEL DO PRIVATE(i,j,k,fx,rhoa,rhob,rs) DEFAULT(NONE)&
     263        12963 : !$OMP SHARED(rho_set,order,e_0,e_ra,e_rb,calc,e_rara,e_rarb,e_rbrb,e_rarara,e_rararb,e_rarbrb,e_rbrbrb)
     264              :       DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
     265              :          DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
     266              :             DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
     267              : 
     268              :                rhoa = rho_set%rhoa(i, j, k)
     269              :                rhob = rho_set%rhob(i, j, k)
     270              :                fx(1) = rhoa + rhob
     271              : 
     272              :                CALL calc_rs(fx(1), rs)
     273              :                CALL calc_fx(rhoa, rhob, fx, ABS(order))
     274              : 
     275              :                IF (calc(0) .AND. calc(1)) THEN
     276              :                   CALL pade_lsd_01(rhoa, rhob, rs, fx, &
     277              :                                    e_0(i, j, k), e_ra(i, j, k), e_rb(i, j, k))
     278              :                ELSE IF (calc(0)) THEN
     279              :                   CALL pade_lsd_0(rhoa, rhob, rs, fx, e_0(i, j, k))
     280              :                ELSE IF (calc(1)) THEN
     281              :                   CALL pade_lsd_1(rhoa, rhob, rs, fx, &
     282              :                                   e_ra(i, j, k), e_rb(i, j, k))
     283              :                END IF
     284              :                IF (calc(2)) THEN
     285              :                   CALL pade_lsd_2(rhoa, rhob, rs, fx, &
     286              :                                   e_rara(i, j, k), e_rarb(i, j, k), e_rbrb(i, j, k))
     287              :                END IF
     288              :                IF (calc(3)) THEN
     289              :                   CALL pade_lsd_3(rhoa, rhob, rs, fx, &
     290              :                                   e_rarara(i, j, k), e_rararb(i, j, k), e_rarbrb(i, j, k), e_rbrbrb(i, j, k))
     291              :                END IF
     292              :             END DO
     293              :          END DO
     294              :       END DO
     295              : 
     296        12963 :    END SUBROUTINE pade_lsd_pw_eval
     297              : 
     298              : ! **************************************************************************************************
     299              : !> \brief ...
     300              : !> \param n ...
     301              : !> \param rho ...
     302              : !> \param rs ...
     303              : !> \param pot ...
     304              : ! **************************************************************************************************
     305         2964 :    SUBROUTINE pade_lda_0(n, rho, rs, pot)
     306              : 
     307              :       INTEGER, INTENT(IN)                                :: n
     308              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, rs
     309              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: pot
     310              : 
     311              :       INTEGER                                            :: ip
     312              :       REAL(KIND=dp)                                      :: epade, p, q
     313              : 
     314              : !$OMP PARALLEL DO PRIVATE(ip,p,q,epade) DEFAULT(NONE)&
     315         2964 : !$OMP SHARED(n,rho,eps_rho,pot,rs)
     316              :       DO ip = 1, n
     317              :          IF (rho(ip) > eps_rho) THEN
     318              :             p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
     319              :             q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
     320              :             epade = -p/q
     321              :             pot(ip) = pot(ip) + epade*rho(ip)
     322              :          END IF
     323              :       END DO
     324              : 
     325         2964 :    END SUBROUTINE pade_lda_0
     326              : 
     327              : ! **************************************************************************************************
     328              : !> \brief ...
     329              : !> \param n ...
     330              : !> \param rho ...
     331              : !> \param rs ...
     332              : !> \param pot ...
     333              : ! **************************************************************************************************
     334            0 :    SUBROUTINE pade_lda_1(n, rho, rs, pot)
     335              : 
     336              :       INTEGER, INTENT(IN)                                :: n
     337              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, rs
     338              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: pot
     339              : 
     340              :       INTEGER                                            :: ip
     341              :       REAL(KIND=dp)                                      :: depade, dpv, dq, epade, p, q
     342              : 
     343              : !$OMP PARALLEL DO PRIVATE(ip,p,q,epade,dpv,dq,depade) DEFAULT(NONE)&
     344            0 : !$OMP SHARED(n,rho,eps_rho,rs,pot)
     345              : 
     346              :       DO ip = 1, n
     347              :          IF (rho(ip) > eps_rho) THEN
     348              : 
     349              :             p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
     350              :             q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
     351              :             epade = -p/q
     352              : 
     353              :             dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs(ip))*rs(ip)
     354              :             dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs(ip))*rs(ip))*rs(ip)
     355              :             depade = f13*rs(ip)*(dpv*q - p*dq)/(q*q)
     356              : 
     357              :             pot(ip) = pot(ip) + epade + depade
     358              : 
     359              :          END IF
     360              :       END DO
     361              : 
     362            0 :    END SUBROUTINE pade_lda_1
     363              : 
     364              : ! **************************************************************************************************
     365              : !> \brief ...
     366              : !> \param n ...
     367              : !> \param rho ...
     368              : !> \param rs ...
     369              : !> \param pot0 ...
     370              : !> \param pot1 ...
     371              : ! **************************************************************************************************
     372        62727 :    SUBROUTINE pade_lda_01(n, rho, rs, pot0, pot1)
     373              : 
     374              :       INTEGER, INTENT(IN)                                :: n
     375              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, rs
     376              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: pot0, pot1
     377              : 
     378              :       INTEGER                                            :: ip
     379              :       REAL(KIND=dp)                                      :: depade, dpv, dq, epade, p, q
     380              : 
     381              : !$OMP PARALLEL DO PRIVATE(ip,p,q,epade,dpv,dq,depade) DEFAULT(NONE)&
     382        62727 : !$OMP SHARED(n,rho,eps_rho,pot0,pot1)
     383              : 
     384              :       DO ip = 1, n
     385              :          IF (rho(ip) > eps_rho) THEN
     386              : 
     387              :             p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
     388              :             q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
     389              :             epade = -p/q
     390              : 
     391              :             dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs(ip))*rs(ip)
     392              :             dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs(ip))*rs(ip))*rs(ip)
     393              :             depade = f13*rs(ip)*(dpv*q - p*dq)/(q*q)
     394              : 
     395              :             pot0(ip) = pot0(ip) + epade*rho(ip)
     396              :             pot1(ip) = pot1(ip) + epade + depade
     397              : 
     398              :          END IF
     399              :       END DO
     400              : 
     401        62727 :    END SUBROUTINE pade_lda_01
     402              : 
     403              : ! **************************************************************************************************
     404              : !> \brief ...
     405              : !> \param n ...
     406              : !> \param rho ...
     407              : !> \param rs ...
     408              : !> \param pot ...
     409              : ! **************************************************************************************************
     410         5606 :    SUBROUTINE pade_lda_2(n, rho, rs, pot)
     411              : 
     412              :       INTEGER, INTENT(IN)                                :: n
     413              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, rs
     414              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: pot
     415              : 
     416              :       INTEGER                                            :: ip
     417              :       REAL(KIND=dp)                                      :: d2p, d2q, dpv, dq, p, q, rsr, t1, t2, t3
     418              : 
     419              : !$OMP PARALLEL DO PRIVATE(ip,p,q,dpv,dq,d2p,d2q,rsr,t1,t2,t3) DEFAULT(NONE)&
     420         5606 : !$OMP SHARED(n,rho,eps_rho,rs)
     421              : 
     422              :       DO ip = 1, n
     423              :          IF (rho(ip) > eps_rho) THEN
     424              : 
     425              :             p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
     426              :             q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
     427              : 
     428              :             dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs(ip))*rs(ip)
     429              :             dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs(ip))*rs(ip))*rs(ip)
     430              : 
     431              :             d2p = 2.0_dp*a2 + 6.0_dp*a3*rs(ip)
     432              :             d2q = 2.0_dp*b2 + (6.0_dp*b3 + 12.0_dp*b4*rs(ip))*rs(ip)
     433              : 
     434              :             rsr = rs(ip)/rho(ip)
     435              :             t1 = (p*dq - dpv*q)/(q*q)
     436              :             t2 = (d2p*q - p*d2q)/(q*q)
     437              :             t3 = (p*dq*dq - dpv*q*dq)/(q*q*q)
     438              : 
     439              :             pot(ip) = pot(ip) - f13*(f23*t1 + f13*t2*rs(ip) + f23*t3*rs(ip))*rsr
     440              : 
     441              :          END IF
     442              :       END DO
     443              : 
     444         5606 :    END SUBROUTINE pade_lda_2
     445              : 
     446              : ! **************************************************************************************************
     447              : !> \brief ...
     448              : !> \param n ...
     449              : !> \param rho ...
     450              : !> \param rs ...
     451              : !> \param pot ...
     452              : ! **************************************************************************************************
     453            0 :    SUBROUTINE pade_lda_3(n, rho, rs, pot)
     454              : 
     455              :       INTEGER, INTENT(IN)                                :: n
     456              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, rs
     457              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: pot
     458              : 
     459              :       INTEGER                                            :: ip
     460              :       REAL(KIND=dp)                                      :: ab1, ab2, ab3, d2p, d2q, d3p, d3q, dpv, &
     461              :                                                             dq, p, q, rsr1, rsr2, rsr3
     462              : 
     463              : !$OMP PARALLEL DO PRIVATE(ip,p,q,dpv,dq,d2p,d2q,d3p,d3q,ab1,ab2,ab3,rsr1,rsr2,rsr3) DEFAULT(NONE)&
     464            0 : !$OMP SHARED(n,rho,eps_rho,rs,pot)
     465              : 
     466              :       DO ip = 1, n
     467              :          IF (rho(ip) > eps_rho) THEN
     468              : 
     469              :             p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
     470              :             q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
     471              : 
     472              :             dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs(ip))*rs(ip)
     473              :             dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs(ip))*rs(ip))*rs(ip)
     474              : 
     475              :             d2p = 2.0_dp*a2 + 6.0_dp*a3*rs(ip)
     476              :             d2q = 2.0_dp*b2 + (6.0_dp*b3 + 12.0_dp*b4*rs(ip))*rs(ip)
     477              : 
     478              :             d3p = 6.0_dp*a3
     479              :             d3q = 6.0_dp*b3 + 24.0_dp*b4*rs(ip)
     480              : 
     481              :             ab1 = (dpv*q - p*dq)/(q*q)
     482              :             ab2 = (d2p*q*q - p*q*d2q - 2.0_dp*dpv*q*dq + 2.0_dp*p*dq*dq)/(q*q*q)
     483              :             ab3 = (d3p*q*q - p*q*d3q - 3.0_dp*dpv*q*d2q + 3.0_dp*p*dq*d2q)/(q*q*q)
     484              :             ab3 = ab3 - 3.0_dp*ab2*dq/q
     485              :             rsr1 = rs(ip)/(rho(ip)*rho(ip))
     486              :             rsr2 = f13*f13*rs(ip)*rsr1
     487              :             rsr3 = f13*rs(ip)*rsr2
     488              :             rsr1 = -f23*f23*f23*rsr1
     489              :             pot(ip) = pot(ip) + rsr1*ab1 + rsr2*ab2 + rsr3*ab3
     490              : 
     491              :          END IF
     492              :       END DO
     493              : 
     494            0 :    END SUBROUTINE pade_lda_3
     495              : 
     496              : ! **************************************************************************************************
     497              : !> \brief ...
     498              : !> \param rhoa ...
     499              : !> \param rhob ...
     500              : !> \param rs ...
     501              : !> \param fx ...
     502              : !> \param pot0 ...
     503              : ! **************************************************************************************************
     504     51414496 :    SUBROUTINE pade_lsd_0(rhoa, rhob, rs, fx, pot0)
     505              : 
     506              :       REAL(KIND=dp), INTENT(IN)                          :: rhoa, rhob, rs
     507              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: fx
     508              :       REAL(KIND=dp), INTENT(INOUT)                       :: pot0
     509              : 
     510              :       REAL(KIND=dp)                                      :: fa0, fa1, fa2, fa3, fb1, fb2, fb3, fb4, &
     511              :                                                             p, q, rhoab
     512              : 
     513     51414496 :       rhoab = rhoa + rhob
     514              : 
     515     51414496 :       IF (rhoab > eps_rho) THEN
     516              : 
     517     48822077 :          fa0 = a0 + fx(1)*da0
     518     48822077 :          fa1 = a1 + fx(1)*da1
     519     48822077 :          fa2 = a2 + fx(1)*da2
     520     48822077 :          fa3 = a3 + fx(1)*da3
     521     48822077 :          fb1 = b1 + fx(1)*db1
     522     48822077 :          fb2 = b2 + fx(1)*db2
     523     48822077 :          fb3 = b3 + fx(1)*db3
     524     48822077 :          fb4 = b4 + fx(1)*db4
     525              : 
     526     48822077 :          p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
     527     48822077 :          q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
     528              : 
     529     48822077 :          pot0 = pot0 - p/q*rhoab
     530              : 
     531              :       END IF
     532              : 
     533     51414496 :    END SUBROUTINE pade_lsd_0
     534              : 
     535              : ! **************************************************************************************************
     536              : !> \brief ...
     537              : !> \param rhoa ...
     538              : !> \param rhob ...
     539              : !> \param rs ...
     540              : !> \param fx ...
     541              : !> \param pota ...
     542              : !> \param potb ...
     543              : ! **************************************************************************************************
     544            0 :    SUBROUTINE pade_lsd_1(rhoa, rhob, rs, fx, pota, potb)
     545              : 
     546              :       REAL(KIND=dp), INTENT(IN)                          :: rhoa, rhob, rs
     547              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: fx
     548              :       REAL(KIND=dp), INTENT(INOUT)                       :: pota, potb
     549              : 
     550              :       REAL(KIND=dp)                                      :: dc, dpv, dq, dr, dx, fa0, fa1, fa2, fa3, &
     551              :                                                             fb1, fb2, fb3, fb4, p, q, rhoab, xp, xq
     552              : 
     553            0 :       rhoab = rhoa + rhob
     554              : 
     555            0 :       IF (rhoab > eps_rho) THEN
     556              : 
     557            0 :          fa0 = a0 + fx(1)*da0
     558            0 :          fa1 = a1 + fx(1)*da1
     559            0 :          fa2 = a2 + fx(1)*da2
     560            0 :          fa3 = a3 + fx(1)*da3
     561            0 :          fb1 = b1 + fx(1)*db1
     562            0 :          fb2 = b2 + fx(1)*db2
     563            0 :          fb3 = b3 + fx(1)*db3
     564            0 :          fb4 = b4 + fx(1)*db4
     565              : 
     566            0 :          p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
     567            0 :          q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
     568            0 :          dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
     569              :          dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
     570            0 :                                    4.0_dp*fb4*rs)*rs)*rs
     571            0 :          xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
     572            0 :          xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
     573              : 
     574            0 :          dr = (dpv*q - p*dq)/(q*q)
     575            0 :          dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx(2)/rhoab
     576            0 :          dc = f13*rs*dr - p/q
     577              : 
     578            0 :          pota = pota + dc - dx*rhob
     579            0 :          potb = potb + dc + dx*rhoa
     580              : 
     581              :       END IF
     582              : 
     583            0 :    END SUBROUTINE pade_lsd_1
     584              : 
     585              : ! **************************************************************************************************
     586              : !> \brief ...
     587              : !> \param rhoa ...
     588              : !> \param rhob ...
     589              : !> \param rs ...
     590              : !> \param fx ...
     591              : !> \param pot0 ...
     592              : !> \param pota ...
     593              : !> \param potb ...
     594              : ! **************************************************************************************************
     595    446926473 :    SUBROUTINE pade_lsd_01(rhoa, rhob, rs, fx, pot0, pota, potb)
     596              : 
     597              :       REAL(KIND=dp), INTENT(IN)                          :: rhoa, rhob, rs
     598              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: fx
     599              :       REAL(KIND=dp), INTENT(INOUT)                       :: pot0, pota, potb
     600              : 
     601              :       REAL(KIND=dp)                                      :: dc, dpv, dq, dr, dx, fa0, fa1, fa2, fa3, &
     602              :                                                             fb1, fb2, fb3, fb4, p, q, rhoab, xp, xq
     603              : 
     604    446926473 :       rhoab = rhoa + rhob
     605              : 
     606    446926473 :       IF (rhoab > eps_rho) THEN
     607              : 
     608    417292816 :          fa0 = a0 + fx(1)*da0
     609    417292816 :          fa1 = a1 + fx(1)*da1
     610    417292816 :          fa2 = a2 + fx(1)*da2
     611    417292816 :          fa3 = a3 + fx(1)*da3
     612    417292816 :          fb1 = b1 + fx(1)*db1
     613    417292816 :          fb2 = b2 + fx(1)*db2
     614    417292816 :          fb3 = b3 + fx(1)*db3
     615    417292816 :          fb4 = b4 + fx(1)*db4
     616              : 
     617    417292816 :          p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
     618    417292816 :          q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
     619    417292816 :          dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
     620              :          dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
     621    417292816 :                                    4.0_dp*fb4*rs)*rs)*rs
     622    417292816 :          xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
     623    417292816 :          xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
     624              : 
     625    417292816 :          dr = (dpv*q - p*dq)/(q*q)
     626    417292816 :          dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx(2)/rhoab
     627    417292816 :          dc = f13*rs*dr - p/q
     628              : 
     629    417292816 :          pot0 = pot0 - p/q*rhoab
     630    417292816 :          pota = pota + dc - dx*rhob
     631    417292816 :          potb = potb + dc + dx*rhoa
     632              : 
     633              :       END IF
     634              : 
     635    446926473 :    END SUBROUTINE pade_lsd_01
     636              : 
     637              : ! **************************************************************************************************
     638              : !> \brief ...
     639              : !> \param rhoa ...
     640              : !> \param rhob ...
     641              : !> \param rs ...
     642              : !> \param fx ...
     643              : !> \param potaa ...
     644              : !> \param potab ...
     645              : !> \param potbb ...
     646              : ! **************************************************************************************************
     647     25165204 :    SUBROUTINE pade_lsd_2(rhoa, rhob, rs, fx, potaa, potab, potbb)
     648              : 
     649              :       REAL(KIND=dp), INTENT(IN)                          :: rhoa, rhob, rs
     650              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: fx
     651              :       REAL(KIND=dp), INTENT(INOUT)                       :: potaa, potab, potbb
     652              : 
     653              :       REAL(KIND=dp)                                      :: d2p, d2q, dpv, dq, dr, drr, dx, dxp, &
     654              :                                                             dxq, dxr, dxx, fa0, fa1, fa2, fa3, &
     655              :                                                             fb1, fb2, fb3, fb4, or, p, q, rhoab, &
     656              :                                                             xp, xq, xt, yt
     657              : 
     658     25165204 :       rhoab = rhoa + rhob
     659              : 
     660     25165204 :       IF (rhoab > eps_rho) THEN
     661              : 
     662     25074450 :          fa0 = a0 + fx(1)*da0
     663     25074450 :          fa1 = a1 + fx(1)*da1
     664     25074450 :          fa2 = a2 + fx(1)*da2
     665     25074450 :          fa3 = a3 + fx(1)*da3
     666     25074450 :          fb1 = b1 + fx(1)*db1
     667     25074450 :          fb2 = b2 + fx(1)*db2
     668     25074450 :          fb3 = b3 + fx(1)*db3
     669     25074450 :          fb4 = b4 + fx(1)*db4
     670              : 
     671     25074450 :          p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
     672     25074450 :          q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
     673              : 
     674     25074450 :          dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
     675              :          dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
     676     25074450 :                                    4.0_dp*fb4*rs)*rs)*rs
     677              : 
     678     25074450 :          d2p = 2.0_dp*fa2 + 6.0_dp*fa3*rs
     679     25074450 :          d2q = 2.0_dp*fb2 + (6.0_dp*fb3 + 12.0_dp*fb4*rs)*rs
     680              : 
     681     25074450 :          xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
     682     25074450 :          xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
     683              : 
     684     25074450 :          dxp = da1 + (2.0_dp*da2 + 3.0_dp*da3*rs)*rs
     685              :          dxq = db1 + (2.0_dp*db2 + (3.0_dp*db3 + &
     686     25074450 :                                     4.0_dp*db4*rs)*rs)*rs
     687              : 
     688     25074450 :          dr = (dpv*q - p*dq)/(q*q)
     689     25074450 :          drr = (d2p*q*q - p*q*d2q - 2.0_dp*dpv*q*dq + 2.0_dp*p*dq*dq)/(q*q*q)
     690     25074450 :          dx = (xp*q - p*xq)/(q*q)
     691     25074450 :          dxx = 2.0_dp*xq*(p*xq - xp*q)/(q*q*q)
     692     25074450 :          dxr = (dxp*q*q + dpv*xq*q - xp*dq*q - p*dxq*q - 2.0_dp*dpv*q*xq + 2.0_dp*p*dq*xq)/(q*q*q)
     693              : 
     694     25074450 :          or = 1.0_dp/rhoab
     695     25074450 :          yt = rhob*or
     696     25074450 :          xt = rhoa*or
     697              : 
     698              :          potaa = potaa + f23*f13*dr*rs*or - f13*f13*drr*rs*rs*or &
     699              :                  + f43*rs*fx(2)*dxr*yt*or &
     700              :                  - 4.0_dp*fx(2)*fx(2)*dxx*yt*yt*or &
     701     25074450 :                  - 4.0_dp*dx*fx(3)*yt*yt*or
     702              :          potab = potab + f23*f13*dr*rs*or - f13*f13*drr*rs*rs*or &
     703              :                  + f23*rs*fx(2)*dxr*(yt - xt)*or &
     704              :                  + 4.0_dp*fx(2)*fx(2)*dxx*xt*yt*or &
     705     25074450 :                  + 4.0_dp*dx*fx(3)*xt*yt*or
     706              :          potbb = potbb + f23*f13*dr*rs*or - f13*f13*drr*rs*rs*or &
     707              :                  - f43*rs*fx(2)*dxr*xt*or &
     708              :                  - 4.0_dp*fx(2)*fx(2)*dxx*xt*xt*or &
     709     25074450 :                  - 4.0_dp*dx*fx(3)*xt*xt*or
     710              : 
     711              :       END IF
     712              : 
     713     25165204 :    END SUBROUTINE pade_lsd_2
     714              : 
     715              : ! **************************************************************************************************
     716              : !> \brief ...
     717              : !> \param rhoa ...
     718              : !> \param rhob ...
     719              : !> \param rs ...
     720              : !> \param fx ...
     721              : !> \param potaaa ...
     722              : !> \param potaab ...
     723              : !> \param potabb ...
     724              : !> \param potbbb ...
     725              : ! **************************************************************************************************
     726            0 :    SUBROUTINE pade_lsd_3(rhoa, rhob, rs, fx, potaaa, potaab, potabb, potbbb)
     727              : 
     728              :       REAL(KIND=dp), INTENT(IN)                          :: rhoa, rhob, rs
     729              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: fx
     730              :       REAL(KIND=dp), INTENT(INOUT)                       :: potaaa, potaab, potabb, potbbb
     731              : 
     732              :       REAL(KIND=dp) :: d2p, d2q, d2xp, d2xq, d3p, d3q, dpv, dq, dr, drr, drrr, dx, dxp, dxq, dxr, &
     733              :          dxrr, dxx, dxxr, dxxx, fa0, fa1, fa2, fa3, fb1, fb2, fb3, fb4, or, p, q, rhoab, xp, xq, &
     734              :          xt, yt
     735              : 
     736            0 :       IF (.NOT. debug_flag) CPABORT("Routine not tested")
     737              : 
     738            0 :       rhoab = rhoa + rhob
     739              : 
     740            0 :       IF (rhoab > eps_rho) THEN
     741              : 
     742            0 :          fa0 = a0 + fx(1)*da0
     743            0 :          fa1 = a1 + fx(1)*da1
     744            0 :          fa2 = a2 + fx(1)*da2
     745            0 :          fa3 = a3 + fx(1)*da3
     746            0 :          fb1 = b1 + fx(1)*db1
     747            0 :          fb2 = b2 + fx(1)*db2
     748            0 :          fb3 = b3 + fx(1)*db3
     749            0 :          fb4 = b4 + fx(1)*db4
     750              : 
     751            0 :          p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
     752            0 :          q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
     753              : 
     754            0 :          dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
     755              :          dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
     756            0 :                                    4.0_dp*fb4*rs)*rs)*rs
     757              : 
     758            0 :          d2p = 2.0_dp*fa2 + 6.0_dp*fa3*rs
     759            0 :          d2q = 2.0_dp*fb2 + (6.0_dp*fb3 + 12.0_dp*fb4*rs)*rs
     760              : 
     761            0 :          d3p = 6.0_dp*fa3
     762            0 :          d3q = 6.0_dp*fb3 + 24.0_dp*fb4*rs
     763              : 
     764            0 :          xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
     765            0 :          xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
     766              : 
     767            0 :          dxp = da1 + (2.0_dp*da2 + 3.0_dp*da3*rs)*rs
     768              :          dxq = db1 + (2.0_dp*db2 + (3.0_dp*db3 + &
     769            0 :                                     4.0_dp*db4*rs)*rs)*rs
     770              : 
     771            0 :          d2xp = 2.0_dp*da2 + 6.0_dp*da3*rs
     772            0 :          d2xq = 2.0_dp*db2 + (6.0_dp*db3 + 12.0_dp*db4*rs)*rs
     773              : 
     774            0 :          dr = (dpv*q - p*dq)/(q*q)
     775            0 :          drr = (d2p*q*q - p*q*d2q - 2.0_dp*dpv*q*dq + 2.0_dp*p*dq*dq)/(q*q*q)
     776              :          drrr = (d3p*q*q*q - 3.0_dp*d2p*dq*q*q + 6.0_dp*dpv*dq*dq*q - 3.0_dp*dpv*d2q*q*q - &
     777            0 :                  6.0_dp*p*dq*dq*dq + 6.0_dp*p*dq*d2q*q - p*d3q*q*q)/(q*q*q*q)
     778            0 :          dx = (xp*q - p*xq)/(q*q)
     779            0 :          dxx = 2.0_dp*xq*(p*xq - xp*q)/(q*q*q)
     780            0 :          dxxx = 6.0_dp*xq*(q*xp*xq - p*xq*xq)/(q*q*q*q)
     781            0 :          dxr = (dxp*q*q + dpv*xq*q - xp*dq*q - p*dxq*q - 2.0_dp*dpv*q*xq + 2.0_dp*p*dq*xq)/(q*q*q)
     782              :          dxxr = 2.0_dp*(2.0_dp*dxq*q*p*xq - dxq*q*q*xp + xq*xq*q*dpv - xq*q*q*dxp + &
     783            0 :                         2.0_dp*xq*q*xp*dq - 3.0_dp*xq*xq*dq*p)/(q*q*q*q)
     784              :          dxrr = (q*q*q*d2xp - 2.0_dp*q*q*dxp*dq - q*q*xp*d2q - q*q*d2p*xq - &
     785              :                  2.0_dp*q*q*dpv*dxq - q*q*p*d2xq + 4.0_dp*dq*q*dpv*xq + 4.0_dp*dq*q*p*dxq + &
     786            0 :                  2.0_dp*dq*dq*q*xp - 6.0_dp*dq*dq*p*xq + 2.0_dp*d2q*q*p*xq)/(q*q*q*q)
     787              : 
     788            0 :          or = 1.0_dp/rhoab
     789            0 :          yt = rhob*or
     790            0 :          xt = rhoa*or
     791              : 
     792              :          potaaa = potaaa + 8.0_dp/27.0_dp*dr*rs*or*or + &
     793              :                   1.0_dp/9.0_dp*drr*rs*rs*or*or + &
     794              :                   1.0_dp/27.0_dp*drrr*rs**3*or*or + &
     795            0 :                   dxr*or*or*yt*rs*(-8.0_dp/3.0_dp*fx(2) + 4.0_dp*fx(3)*yt)
     796            0 :          potaab = potaab + 0.0_dp
     797            0 :          potabb = potabb + 0.0_dp
     798            0 :          potbbb = potbbb + 0.0_dp
     799              : 
     800              :       END IF
     801              : 
     802            0 :    END SUBROUTINE pade_lsd_3
     803              : 
     804              : ! **************************************************************************************************
     805              : !> \brief ...
     806              : !> \param rho_a ...
     807              : !> \param rho_b ...
     808              : !> \param fxc_aa ...
     809              : !> \param fxc_ab ...
     810              : !> \param fxc_bb ...
     811              : ! **************************************************************************************************
     812            2 :    SUBROUTINE pade_fxc_eval(rho_a, rho_b, fxc_aa, fxc_ab, fxc_bb)
     813              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: rho_a, rho_b
     814              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: fxc_aa, fxc_ab, fxc_bb
     815              : 
     816              :       INTEGER                                            :: i, j, k
     817              :       INTEGER, DIMENSION(2, 3)                           :: bo
     818              :       REAL(KIND=dp)                                      :: eaa, eab, ebb, rhoa, rhob, rs
     819              :       REAL(KIND=dp), DIMENSION(4)                        :: fx
     820              : 
     821           20 :       bo(1:2, 1:3) = rho_a%pw_grid%bounds_local(1:2, 1:3)
     822              : !$OMP PARALLEL DO PRIVATE(i,j,k,fx,rhoa,rhob,rs,eaa,eab,ebb) DEFAULT(NONE)&
     823            2 : !$OMP SHARED(bo,rho_a,rho_b,fxc_aa,fxc_ab,fxc_bb)
     824              :       DO k = bo(1, 3), bo(2, 3)
     825              :          DO j = bo(1, 2), bo(2, 2)
     826              :             DO i = bo(1, 1), bo(2, 1)
     827              : 
     828              :                rhoa = rho_a%array(i, j, k)
     829              :                rhob = rho_b%array(i, j, k)
     830              :                fx(1) = rhoa + rhob
     831              : 
     832              :                CALL calc_rs(fx(1), rs)
     833              :                CALL calc_fx(rhoa, rhob, fx, 2)
     834              : 
     835              :                eaa = 0.0_dp; eab = 0.0_dp; ebb = 0.0_dp
     836              :                CALL pade_lsd_2(rhoa, rhob, rs, fx, eaa, eab, ebb)
     837              : 
     838              :                fxc_aa%array(i, j, k) = eaa
     839              :                fxc_ab%array(i, j, k) = eab
     840              :                fxc_bb%array(i, j, k) = ebb
     841              : 
     842              :             END DO
     843              :          END DO
     844              :       END DO
     845              : 
     846            2 :    END SUBROUTINE pade_fxc_eval
     847              : 
     848              : END MODULE xc_pade
     849              : 
        

Generated by: LCOV version 2.0-1