LCOV - code coverage report
Current view: top level - src/xc - xc_perdew86.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 69.2 % 78 54
Test Date: 2025-07-25 12:55:17 Functions: 71.4 % 7 5

            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 Perdew Correlation from 1986
      10              : !> \par History
      11              : !>      JGH (26.02.2003) : OpenMP enabled
      12              : !>      fawzi (04.2004)  : adapted to the new xc interface
      13              : !> \author JGH (03.03.2002)
      14              : ! **************************************************************************************************
      15              : MODULE xc_perdew86
      16              : 
      17              :    USE input_section_types,             ONLY: section_vals_type
      18              :    USE kinds,                           ONLY: dp
      19              :    USE xc_derivative_desc,              ONLY: deriv_norm_drho,&
      20              :                                               deriv_rho
      21              :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      22              :                                               xc_dset_get_derivative
      23              :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      24              :                                               xc_derivative_type
      25              :    USE xc_functionals_utilities,        ONLY: calc_rs_pw,&
      26              :                                               set_util
      27              :    USE xc_input_constants,              ONLY: pz_orig
      28              :    USE xc_perdew_zunger,                ONLY: pz_lda_eval
      29              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      30              :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      31              :                                               xc_rho_set_type
      32              : #include "../base/base_uses.f90"
      33              : 
      34              :    IMPLICIT NONE
      35              : 
      36              :    PRIVATE
      37              : 
      38              : ! *** Global parameters ***
      39              : 
      40              :    REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp
      41              :    REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
      42              :                                f23 = 2.0_dp*f13, &
      43              :                                f43 = 4.0_dp*f13, &
      44              :                                f53 = 5.0_dp*f13, &
      45              :                                f76 = 7.0_dp/6.0_dp, &
      46              :                                frs = 1.6119919540164696407_dp, &
      47              :                                fpe = 0.19199566167376364_dp
      48              : 
      49              :    PUBLIC :: p86_lda_info, p86_lda_eval
      50              : 
      51              :    REAL(KIND=dp) :: eps_rho
      52              :    LOGICAL :: debug_flag
      53              : 
      54              :    REAL(KIND=dp), PARAMETER :: a = 0.023266_dp, &
      55              :                                b = 7.389e-6_dp, &
      56              :                                c = 8.723_dp, &
      57              :                                d = 0.472_dp, &
      58              :                                pc1 = 0.001667_dp, &
      59              :                                pc2 = 0.002568_dp, &
      60              :                                pci = pc1 + pc2
      61              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_perdew86'
      62              : 
      63              : CONTAINS
      64              : 
      65              : ! **************************************************************************************************
      66              : !> \brief ...
      67              : !> \param cutoff ...
      68              : !> \param debug ...
      69              : ! **************************************************************************************************
      70            4 :    SUBROUTINE p86_init(cutoff, debug)
      71              : 
      72              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
      73              :       LOGICAL, INTENT(IN), OPTIONAL                      :: debug
      74              : 
      75            4 :       eps_rho = cutoff
      76            4 :       CALL set_util(cutoff)
      77              : 
      78            4 :       IF (PRESENT(debug)) THEN
      79            0 :          debug_flag = debug
      80              :       ELSE
      81            4 :          debug_flag = .FALSE.
      82              :       END IF
      83              : 
      84            4 :    END SUBROUTINE p86_init
      85              : 
      86              : ! **************************************************************************************************
      87              : !> \brief ...
      88              : !> \param reference ...
      89              : !> \param shortform ...
      90              : !> \param needs ...
      91              : !> \param max_deriv ...
      92              : ! **************************************************************************************************
      93            9 :    SUBROUTINE p86_lda_info(reference, shortform, needs, max_deriv)
      94              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      95              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      96              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      97              : 
      98            9 :       IF (PRESENT(reference)) THEN
      99            1 :          reference = "J. P. Perdew, Phys. Rev. B, 33, 8822 (1986) {LDA version}"
     100              :       END IF
     101            9 :       IF (PRESENT(shortform)) THEN
     102            1 :          shortform = "Perdew 1986 correlation energy functional {LDA}"
     103              :       END IF
     104            9 :       IF (PRESENT(needs)) THEN
     105            8 :          needs%rho = .TRUE.
     106            8 :          needs%norm_drho = .TRUE.
     107              :       END IF
     108            9 :       IF (PRESENT(max_deriv)) max_deriv = 3
     109              : 
     110            9 :    END SUBROUTINE p86_lda_info
     111              : 
     112              : ! **************************************************************************************************
     113              : !> \brief ...
     114              : !> \param rho_set ...
     115              : !> \param deriv_set ...
     116              : !> \param order ...
     117              : !> \param p86_params ...
     118              : ! **************************************************************************************************
     119            4 :    SUBROUTINE p86_lda_eval(rho_set, deriv_set, order, p86_params)
     120              : 
     121              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     122              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     123              :       INTEGER, INTENT(IN)                                :: order
     124              :       TYPE(section_vals_type), POINTER                   :: p86_params
     125              : 
     126              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'p86_lda_eval'
     127              : 
     128              :       INTEGER                                            :: handle, m, npoints
     129              :       INTEGER, DIMENSION(2, 3)                           :: bo
     130              :       REAL(KIND=dp)                                      :: drho_cutoff, rho_cutoff
     131            4 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: rs
     132            4 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
     133            4 :          e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
     134            4 :          e_rho_rho_rho, grho, rho
     135              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     136              : 
     137            4 :       CALL timeset(routineN, handle)
     138            4 :       NULLIFY (rho, e_0, e_rho, e_ndrho, &
     139            4 :                e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
     140            4 :                e_rho_rho_rho, e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho)
     141              : 
     142              :       ! calculate the perdew_zunger correlation
     143            4 :       CALL pz_lda_eval(pz_orig, rho_set, deriv_set, order, p86_params)
     144              : 
     145              :       CALL xc_rho_set_get(rho_set, rho=rho, &
     146              :                           norm_drho=grho, local_bounds=bo, rho_cutoff=rho_cutoff, &
     147            4 :                           drho_cutoff=drho_cutoff)
     148            4 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     149            4 :       CALL p86_init(rho_cutoff)
     150            4 :       m = ABS(order)
     151              : 
     152           12 :       ALLOCATE (rs(npoints))
     153              : 
     154            4 :       CALL calc_rs_pw(rho, rs, npoints)
     155            4 :       IF (order >= 0) THEN
     156              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     157            4 :                                          allocate_deriv=.TRUE.)
     158            4 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     159              : 
     160            4 :          CALL p86_u_0(rho, rs, grho, e_0, npoints)
     161              :       END IF
     162            4 :       IF (order >= 1 .OR. order == -1) THEN
     163              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     164            4 :                                          allocate_deriv=.TRUE.)
     165            4 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     166              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     167            4 :                                          allocate_deriv=.TRUE.)
     168            4 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     169              : 
     170              :          CALL p86_u_1(rho, grho, rs, e_rho, &
     171            4 :                       e_ndrho, npoints)
     172              :       END IF
     173            4 :       IF (order >= 2 .OR. order == -2) THEN
     174              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     175            0 :                                          allocate_deriv=.TRUE.)
     176            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     177              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
     178            0 :                                          allocate_deriv=.TRUE.)
     179            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
     180              :          deriv => xc_dset_get_derivative(deriv_set, &
     181            0 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     182            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     183              : 
     184              :          CALL p86_u_2(rho, grho, rs, e_rho_rho, &
     185            0 :                       e_rho_ndrho, e_ndrho_ndrho, npoints)
     186              :       END IF
     187            4 :       IF (order >= 3 .OR. order == -3) THEN
     188              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
     189            0 :                                          allocate_deriv=.TRUE.)
     190            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     191              :          deriv => xc_dset_get_derivative(deriv_set, &
     192            0 :                                          [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.TRUE.)
     193            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
     194              :          deriv => xc_dset_get_derivative(deriv_set, &
     195            0 :                                          [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     196            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
     197              :          deriv => xc_dset_get_derivative(deriv_set, &
     198            0 :                                          [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     199            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
     200              : 
     201              :          CALL p86_u_3(rho, grho, rs, e_rho_rho_rho, &
     202              :                       e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
     203            0 :                       npoints)
     204              :       END IF
     205            4 :       IF (order > 3 .OR. order < -3) THEN
     206            0 :          CPABORT("derivatives bigger than 3 not implemented")
     207              :       END IF
     208            4 :       DEALLOCATE (rs)
     209            4 :       CALL timestop(handle)
     210              : 
     211            8 :    END SUBROUTINE p86_lda_eval
     212              : 
     213              : ! **************************************************************************************************
     214              : !> \brief ...
     215              : !> \param rho ...
     216              : !> \param rs ...
     217              : !> \param grho ...
     218              : !> \param e_0 ...
     219              : !> \param npoints ...
     220              : ! **************************************************************************************************
     221            4 :    SUBROUTINE p86_u_0(rho, rs, grho, e_0, npoints)
     222              : 
     223              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, rs, grho
     224              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_0
     225              :       INTEGER, INTENT(in)                                :: npoints
     226              : 
     227              :       INTEGER                                            :: ip
     228              :       REAL(KIND=dp)                                      :: cr, ep, g, or, phi, r, x
     229              : 
     230              : !$OMP PARALLEL DO PRIVATE(ip,g,r,x,or,cr,phi,ep) DEFAULT(NONE)&
     231            4 : !$OMP SHARED(npoints,rho,eps_rho,grho,rs,e_0)
     232              :       DO ip = 1, npoints
     233              : 
     234              :          IF (rho(ip) > eps_rho) THEN
     235              :             g = grho(ip)
     236              :             r = rs(ip)
     237              :             x = r*frs
     238              :             or = 1.0_dp/rho(ip)
     239              :             cr = pc1 + (pc2 + a*r + b*r*r)/(1.0_dp + c*r + d*r*r + 1.e4_dp*b*r*r*r)
     240              :             phi = fpe*pci/cr*g*SQRT(x)*or
     241              :             ep = EXP(-phi)
     242              :             e_0(ip) = e_0(ip) + x*or*g*g*cr*ep
     243              :          END IF
     244              : 
     245              :       END DO
     246              : 
     247            4 :    END SUBROUTINE p86_u_0
     248              : 
     249              : ! **************************************************************************************************
     250              : !> \brief ...
     251              : !> \param rho ...
     252              : !> \param grho ...
     253              : !> \param rs ...
     254              : !> \param e_rho ...
     255              : !> \param e_ndrho ...
     256              : !> \param npoints ...
     257              : ! **************************************************************************************************
     258            4 :    SUBROUTINE p86_u_1(rho, grho, rs, e_rho, e_ndrho, npoints)
     259              : 
     260              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, grho, rs
     261              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho, e_ndrho
     262              :       INTEGER, INTENT(in)                                :: npoints
     263              : 
     264              :       INTEGER                                            :: ip
     265              :       REAL(KIND=dp)                                      :: cr, dcr, dphig, dphir, dpv, dq, ep, ff, &
     266              :                                                             g, or, p, phi, q, r, x
     267              : 
     268              : !$OMP PARALLEL DO PRIVATE(ip,g,r,x,or,p,dpv,q,dq,cr,dcr,dphig,phi,dphir,ep,ff) DEFAULT(NONE)&
     269            4 : !$OMP SHARED(npoints,rho,eps_rho,grho,rs,e_rho,e_ndrho)
     270              :       DO ip = 1, npoints
     271              : 
     272              :          IF (rho(ip) > eps_rho) THEN
     273              :             g = grho(ip)
     274              :             r = rs(ip)
     275              :             x = r*frs
     276              :             or = 1.0_dp/rho(ip)
     277              :             p = pc2 + a*r + b*r*r
     278              :             dpv = a + 2.0_dp*b*r
     279              :             q = 1.0_dp + c*r + d*r*r + 1.e4_dp*b*r*r*r
     280              :             dq = c + 2.0_dp*d*r + 3.e4_dp*b*r*r
     281              :             cr = pc1 + p/q
     282              :             dcr = (dpv*q - p*dq)/(q*q)*(-f13*r*or)
     283              :             dphig = fpe*pci/cr*SQRT(x)*or
     284              :             phi = dphig*g
     285              :             dphir = -phi*(dcr/cr + f76*or)
     286              :             ep = EXP(-phi)
     287              :             ff = x*or*g*ep
     288              :             e_rho(ip) = e_rho(ip) + ff*g*dcr - ff*g*cr*dphir - ff*g*cr*f43*or
     289              :             e_ndrho(ip) = e_ndrho(ip) + ff*cr*(2.0_dp - g*dphig)
     290              :          END IF
     291              : 
     292              :       END DO
     293              : 
     294            4 :    END SUBROUTINE p86_u_1
     295              : 
     296              : ! **************************************************************************************************
     297              : !> \brief ...
     298              : !> \param rho ...
     299              : !> \param grho ...
     300              : !> \param rs ...
     301              : !> \param e_rho_rho ...
     302              : !> \param e_rho_ndrho ...
     303              : !> \param e_ndrho_ndrho ...
     304              : !> \param npoints ...
     305              : ! **************************************************************************************************
     306            0 :    SUBROUTINE p86_u_2(rho, grho, rs, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
     307              :                       npoints)
     308              : 
     309              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, grho, rs
     310              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
     311              :       INTEGER, INTENT(in)                                :: npoints
     312              : 
     313              :       INTEGER                                            :: ip
     314              :       REAL(KIND=dp)                                      :: cr, d2cr, d2p, d2phir, d2q, dcr, dphig, &
     315              :                                                             dphigr, dphir, dpv, dq, ep, g, or, p, &
     316              :                                                             phi, q, r, x
     317              : 
     318              : !$OMP PARALLEL DO PRIVATE(ip,x,r,cr,phi,ep,g,or,p,q,dpv,dq,dphir,dcr) &
     319              : !$OMP             PRIVATE(dphig,dphigr,d2phir,d2cr,d2p,d2q) DEFAULT(NONE) &
     320            0 : !$OMP SHARED(npoints,rho,eps_rho,grho,rs,e_rho_rho,e_rho_ndrho,e_ndrho_ndrho)
     321              :       DO ip = 1, npoints
     322              : 
     323              :          IF (rho(ip) > eps_rho) THEN
     324              :             g = grho(ip)
     325              :             r = rs(ip)
     326              :             x = r*frs
     327              :             or = 1.0_dp/rho(ip)
     328              :             p = pc2 + a*r + b*r*r
     329              :             dpv = a + 2.0_dp*b*r
     330              :             d2p = 2.0_dp*b
     331              :             q = 1.0_dp + c*r + d*r*r + 1.e4_dp*b*r*r*r
     332              :             dq = c + 2.0_dp*d*r + 3.e4_dp*b*r*r
     333              :             d2q = 2.0_dp*d + 6.e4_dp*b*r
     334              :             cr = pc1 + p/q
     335              :             dcr = (dpv*q - p*dq)/(q*q)*(-f13*r*or)
     336              :             d2cr = (d2p*q*q - p*q*d2q - 2*dpv*dq*q + 2*p*dq*dq)/(q*q*q)*(f13*r*or)**2 + &
     337              :                    (dpv*q - p*dq)/(q*q)*f13*f43*r*or*or
     338              :             dphig = fpe*pci/cr*SQRT(x)*or
     339              :             phi = dphig*g
     340              :             dphir = -phi*(dcr/cr + f76*or)
     341              :             d2phir = -dphir*(dcr/cr + f76*or) - &
     342              :                      phi*((d2cr*cr - dcr*dcr)/(cr*cr) - f76*or*or)
     343              :             dphigr = -dphig*(dcr/cr + f76*or)
     344              :             ep = EXP(-phi)
     345              :             e_rho_rho(ip) = e_rho_rho(ip) + x*or*ep*g*g* &
     346              :                             (-f43*or*dcr + d2cr - dcr*dphir + &
     347              :                              f43*or*cr*dphir - dcr*dphir - cr*d2phir + cr*dphir*dphir + &
     348              :                              f43*or*(7.*f13*or*cr - dcr + cr*dphir))
     349              :             e_rho_ndrho(ip) = e_rho_ndrho(ip) + x*or*ep*g* &
     350              :                               (-2*f43*cr*or + 2*dcr - 2*cr*dphir + f43*or*g*cr*dphig - &
     351              :                                g*dcr*dphig + g*cr*dphir*dphig - g*cr*dphigr)
     352              :             e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + x*or*ep*cr* &
     353              :                                 (2.0_dp - 4.0_dp*g*dphig + g*g*dphig*dphig)
     354              :          END IF
     355              : 
     356              :       END DO
     357              : 
     358            0 :    END SUBROUTINE p86_u_2
     359              : 
     360              : ! **************************************************************************************************
     361              : !> \brief ...
     362              : !> \param rho ...
     363              : !> \param grho ...
     364              : !> \param rs ...
     365              : !> \param e_rho_rho_rho ...
     366              : !> \param e_rho_rho_ndrho ...
     367              : !> \param e_rho_ndrho_ndrho ...
     368              : !> \param e_ndrho_ndrho_ndrho ...
     369              : !> \param npoints ...
     370              : ! **************************************************************************************************
     371            0 :    SUBROUTINE p86_u_3(rho, grho, rs, e_rho_rho_rho, &
     372              :                       e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
     373              :                       npoints)
     374              : 
     375              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, grho, rs
     376              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho_rho, e_rho_rho_ndrho, &
     377              :                                                             e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
     378              :       INTEGER, INTENT(in)                                :: npoints
     379              : 
     380              :       INTEGER                                            :: ip
     381              :       REAL(KIND=dp) :: cr, d2cr, d2p, d2phir, d2phirg, d2pq, d2q, d2z, d3cr, d3phir, d3pq, d3q, &
     382              :          d3z, dcr, dphig, dphigr, dphir, dpq, dpv, dq, dz, ep, g, or, oz, p, phi, pq, q, r, x
     383              : 
     384              : !$OMP PARALLEL DO PRIVATE(ip,x, r, cr, phi, ep, g, or, p, q, dpv, dq, dphir, dcr, dphig) &
     385              : !$OMP             PRIVATE(dphigr, d2phir, d3phir, d2cr, d3cr, d2p, d2q, d2phirg, d3q) &
     386              : !$OMP             PRIVATE(pq, dpq, d2pq, d3pq, oz, dz, d2z, d3z) DEFAULT(NONE) &
     387            0 : !$OMP             SHARED(npoints,rho,eps_rho,grho,e_rho_rho_rho,e_rho_rho_ndrho,e_rho_ndrho_ndrho,e_ndrho_ndrho_ndrho)
     388              :       DO ip = 1, npoints
     389              : 
     390              :          IF (rho(ip) > eps_rho) THEN
     391              :             g = grho(ip)
     392              :             r = rs(ip)
     393              :             x = r*frs
     394              :             or = 1.0_dp/rho(ip)
     395              :             p = pc2 + a*r + b*r*r
     396              :             dpv = a + 2.0_dp*b*r
     397              :             d2p = 2.0_dp*b
     398              :             q = 1.0_dp + c*r + d*r*r + 1.e4_dp*b*r*r*r
     399              :             dq = c + 2.0_dp*d*r + 3.e4_dp*b*r*r
     400              :             d2q = 2.0_dp*d + 6.e4*b*r
     401              :             d3q = 6.e4*b
     402              :             pq = p/q
     403              :             dpq = (dpv*q - p*dq)/(q*q)
     404              :             d2pq = (d2p*q*q - 2*dpv*dq*q + 2*p*dq*dq - p*d2q*q)/(q*q*q)
     405              :             d3pq = -(3*d2p*dq*q*q - 6*dpv*dq*dq*q + 3*dpv*d2q*q*q + 6*p*dq*dq*dq - 6*p*dq*d2q*q &
     406              :                      + p*d3q*q*q)/(q*q*q*q)
     407              :             cr = pc1 + pq
     408              :             dcr = dpq*(-f13*r*or)
     409              :             d2cr = d2pq*f13*f13*r*r*or*or + dpq*f13*f43*r*or*or
     410              :             d3cr = d3pq*(-f13*r*or)**3 + 3*d2pq*(-f13*f13*f43*r*r*or*or*or) + &
     411              :                    dpq*(-f13*f43*f13*7*r*or*or*or)
     412              :             oz = SQRT(x)*or/cr
     413              :             dz = dcr/cr + f76*or
     414              :             d2z = d2cr/cr + 2*f76*dcr/cr*or + f76/6.*or*or
     415              :             d3z = d3cr/cr + 3*f76*d2cr/cr*or + 3*f76/6.*dcr/cr*or*or - 5*f76/36.*or*or*or
     416              :             dphig = fpe*pci*oz
     417              :             phi = dphig*g
     418              :             dphir = -phi*dz
     419              :             dphigr = -dphig*dz
     420              :             d2phir = -phi*(d2z - 2*dz*dz)
     421              :             d3phir = -phi*(d3z - 6*d2z*dz + 6*dz*dz*dz)
     422              :             d2phirg = -dphigr*dz - &
     423              :                       dphig*((d2cr*cr - dcr*dcr)/(cr*cr) - f76*or*or)
     424              :             ep = EXP(-phi)
     425              :             e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
     426              :                                 + g*g*x*or*ep*(-280./27.*or*or*or*cr + 3*28./9.*or*or*dcr + &
     427              :                                                3*28./9.*or*or*cr*(-dphir) - 4*or*d2cr - 8*or*dcr*(-dphir) - &
     428              :                                                4*or*cr*(-d2phir + dphir*dphir) + d3cr + 3*d2cr*(-dphir) + &
     429              :                                                3*dcr*(-d2phir + dphir*dphir) + cr*(-d3phir + 3*dphir*d2phir - &
     430              :                                                                                    dphir**3))
     431              :             e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
     432              :                                   + 2.*x*or*ep*g*(-f43*or*dcr + d2cr - dcr*dphir + &
     433              :                                                   f43*or*cr*dphir - dcr*dphir - cr*d2phir + cr*dphir*dphir + &
     434              :                                                   f43*or*(7.*f13*or*cr - dcr + cr*dphir)) - &
     435              :                                   dphig*x*or*ep*g*g*(-f43*or*dcr + d2cr - dcr*dphir + &
     436              :                                                      f43*or*cr*dphir - dcr*dphir - cr*d2phir + cr*dphir*dphir + &
     437              :                                                      f43*or*(7.*f13*or*cr - dcr + cr*dphir)) + &
     438              :                                   x*or*ep*g*g*(-dcr*dphigr + f43*or*cr*dphigr - dcr*dphigr - cr*d2phirg + &
     439              :                                                2.*cr*dphigr*dphir + f43*or*cr*dphigr)
     440              :             e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
     441              :                                     + x*or*ep*(-2*f43*cr*or + 2*dcr - 2*cr*dphir + f43*or*g*cr*dphig - &
     442              :                                                g*dcr*dphig + g*cr*dphir*dphig - g*cr*dphigr) + &
     443              :                                     x*or*ep*g*(-2*cr*dphigr + f43*or*cr*dphig - &
     444              :                                                dcr*dphig + cr*dphir*dphig + g*cr*dphigr*dphig - cr*dphigr) - &
     445              :                                     x*or*ep*g*dphig*(-2*f43*cr*or + 2*dcr - 2*cr*dphir + f43*or*g*cr*dphig - &
     446              :                                                      g*dcr*dphig + g*cr*dphir*dphig - g*cr*dphigr)
     447              :             e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) &
     448              :                                       + x*or*ep*cr*dphig*(-6.0_dp + 6.0_dp*g*dphig - g*g*dphig*dphig)
     449              :          END IF
     450              : 
     451              :       END DO
     452              : 
     453            0 :    END SUBROUTINE p86_u_3
     454              : 
     455              : END MODULE xc_perdew86
     456              : 
        

Generated by: LCOV version 2.0-1