Question: libxc integration in two-center code

[vanderhe]

(this is more of a question)

For fun I replaced the current LDA and PBE routines of sktwocnt by their libxc counterpart, resulting in:

  use, intrinsic :: iso_c_binding, only : c_size_t
  use xc_f90_lib_m, only : xc_f90_func_t, xc_f90_func_info_t, xc_f90_func_init,&
      & xc_f90_func_get_info, xc_f90_lda_vxc, xc_f90_gga_vxc, xc_f90_func_end, XC_LDA_X,&
      & XC_LDA_C_PW, XC_GGA_X_PBE, XC_GGA_C_PBE, XC_UNPOLARIZED

  subroutine getxcpot_ldapw91(rho4pi, xcpot)
    real(dp), intent(in) :: rho4pi(:)
    real(dp), intent(out) :: xcpot(:)

    !> libxc related objects
    type(xc_f90_func_t) :: xcfunc_x, xcfunc_c
    type(xc_f90_func_info_t) :: xcinfo

    !> density with libxc compatible normalization
    real(dp), allocatable :: rho(:)

    !> exchange and correlation potential on grid
    real(dp), allocatable :: vx(:), vc(:)

    !> number of density grid points
    integer(c_size_t) :: nn

    call xc_f90_func_init(xcfunc_x, XC_LDA_X, XC_UNPOLARIZED)
    xcinfo = xc_f90_func_get_info(xcfunc_x)
    call xc_f90_func_init(xcfunc_c, XC_LDA_C_PW, XC_UNPOLARIZED)
    xcinfo = xc_f90_func_get_info(xcfunc_x)

    nn = size(rho4pi)
    allocate(vx(nn), vc(nn))

    rho = rho4pi * rec4pi

    call xc_f90_lda_vxc(xcfunc_x, nn, rho, vx)
    call xc_f90_lda_vxc(xcfunc_c, nn, rho, vc)

    xcpot(:) = vx + vc

    call xc_f90_func_end(xcfunc_x)
    call xc_f90_func_end(xcfunc_c)

  end subroutine getxcpot_ldapw91

  subroutine getxcpot_ggapbe(rho4pi, absgr4pi, laplace4pi, gr_grabsgr4pi, xcpot)
    real(dp), intent(in) :: rho4pi(:)
    real(dp), intent(in) :: absgr4pi(:), laplace4pi(:), gr_grabsgr4pi(:)
    real(dp), intent(out) :: xcpot(:)

    !> libxc related objects
    type(xc_f90_func_t) :: xcfunc_x, xcfunc_c
    type(xc_f90_func_info_t) :: xcinfo

    !> density with libxc compatible normalization
    real(dp), allocatable :: rho(:)

    !> contracted gradients of the density
    real(dp), allocatable :: sigma(:)

    !> exchange and correlation potential on grid
    real(dp), allocatable :: vx(:), vc(:)

    !> first partial derivative of the energy per unit volume in terms of sigma
    real(dp), allocatable :: vxsigma(:), vcsigma(:)

    !> number of density grid points
    integer(c_size_t) :: nn

    nn = size(rho4pi)
    allocate(vx(nn), vc(nn), vxsigma(nn), vcsigma(nn))

    rho = rho4pi * rec4pi
    sigma = (absgr4pi * rec4pi)**2

    call xc_f90_func_init(xcfunc_x, XC_GGA_X_PBE, XC_UNPOLARIZED)
    xcinfo = xc_f90_func_get_info(xcfunc_x)
    call xc_f90_func_init(xcfunc_c, XC_GGA_C_PBE, XC_UNPOLARIZED)
    xcinfo = xc_f90_func_get_info(xcfunc_x)

    call xc_f90_gga_vxc(xcfunc_x, nn, rho, sigma, vx, vxsigma)
    call xc_f90_gga_vxc(xcfunc_c, nn, rho, sigma, vc, vcsigma)

    xcpot(:) = vx + vc

    call xc_f90_func_end(xcfunc_x)
    call xc_f90_func_end(xcfunc_c)

  end subroutine getxcpot_ggapbe

The values for vx(:) and vc(:) however are very different (resulting in deviations 1e-7 a.u. for LDA and 1e-2 a.u. for PBE Hamiltonian matrix elements). Maybe I screwed something up with the normalization of the density (-gradient) or is this a known issue due to possibly different implementations (Burke's Fortran implementation <--> libxc)?