subroutine corggapbe2 (lpot, np, nilrho, rs, zet, g91, t, uu, vv, ww, &
                       ec, ecrs, eczet, h, dvcup, dvcdn, bet)

   use Vartypes
   implicit real*8 (a-h,o-z)

   integer(KINT), intent(IN) :: np
   real(KREAL), intent(IN)   :: rs(np), zet(np), t(np), uu(np), vv(np), ww(np)
   real(KREAL), intent(IN)   :: ec(np), ecrs(np), eczet(np), g91(np), bet
   real(KREAL), intent(OUT)  :: h(np), dvcup(np), dvcdn(np)
   logical, intent(IN)       :: lpot, nilrho(np)

!----------------------------------------------------------------------
!   adaptation of official PBE correlation code for use within ADF:
!   - simplification of t-dependence from:
!       1 + delta*t2*[1+A*t2]/[1+A*t2+A2*t4]
!     to:
!       1 + delta*t2*1/[1+A*t2]
!     analogous to exchange-part
!   adaptation: M. Swart, 2008
!
!   input: rs=seitz radius=(3/4pi rho)^(1/3)
!        : zet=relative spin polarization = (rhoup-rhodn)/rho
!        : t=abs(grad rho)/(rho*2.*ks*g)  -- only needed for pbe
!        : uu=(grad rho)*grad(abs(grad rho))/(rho**2 * (2*ks*g)**3)
!        : vv=(laplacian rho)/(rho * (2*ks*g)**2)
!        : ww=(grad rho)*(grad zet)/(rho * (2*ks*g)**2
!        :  uu,vv,ww, only needed for pbe potential
!        : lpot=flag to do potential
!   output: ec=lsd correlation energy from [a]
!         : vcup=lsd up correlation potential
!         : vcdn=lsd dn correlation potential
!         : h=nonlocal part of correlation energy per electron
!         : dvcup=nonlocal correction to vcup
!         : dvcdn=nonlocal correction to vcdn
!
!  references:
!  [a] j.p. perdew, k. burke, and m. ernzerhof,
!      phys. rev. lett. {\bf 77}, 3865 (1996).
!  [b] j. p. perdew, k. burke, and y. wang,
!      phys. rev. b {\bf 54}, 16533 (1996).
!  [c] j. p. perdew and y. wang, phys. rev. b {\bf 45}, 13244 (1992).
!  [d] m. swart, m. sola, f.m. bickelhaupt, j. chem. phys.
!      {\bf 131}, 094103 (2009).
!
!  original code: K. Burke, July 1996
!                 (cleaned by tof90 utility)
!  adapted for ADF: M. Swart, November 2007
!----------------------------------------------------------------------

!  -------------------------------------------------------------------
!  numbers for construction of pbe
!       gamma=(1-log(2))/pi^2
!       bet=coefficient in gradient expansion for correlation, [a](4).
!       eta=small number to stop d phi/ dzeta from blowing up at
!           |zeta|=1.
!  -------------------------------------------------------------------

   real(KREAL), parameter :: ZERO  = 0.0_KREAL
   real(KREAL), parameter :: ONE   = 1.0_KREAL
   real(KREAL), parameter :: TWO   = 2.0_KREAL
   real(KREAL), parameter :: THREE = 3.0_KREAL
   real(KREAL), parameter :: FOUR  = 4.0_KREAL
   real(KREAL), parameter :: SIX   = 6.0_KREAL
   real(KREAL), parameter :: SEVEN = 7.0_KREAL

   real(KREAL), parameter :: THRD2  =  TWO / THREE
   real(KREAL), parameter :: SIXTHM = -ONE / SIX
   real(KREAL), parameter :: gamma  =  0.0310906908696548950349408637127_KREAL
   real(KREAL), parameter :: eta    =  1.0e-12_KREAL

!  -------------------------------
!  pbe correlation energy
!  g=phi(zeta), given after [a](3)
!  delt=bet/gamma
!  b=a of [a](8)
!  -------------------------------

   integer(KINT) :: i
   real(KREAL)   :: delt
   delt = bet/gamma

   lpi: &         
   do i=1,np

     if (nilrho(i)) then
        h(i) = zero
        if (lpot) then
          dvcup(i) = zero
          dvcdn(i) = zero
        end if
        cycle lpi 
     end if 

!    g   = ((ONE + zet(i))**THRD2 + (ONE - zet(i))**THRD2)/TWO
!    the g has been calculated already before and put into g91-array

     g3  = g91(i)**THREE
     pon = -ec(i)/(g3*gamma)
     b   = delt/(exp(pon) - ONE)
     t2  = t(i)*t(i)
     q2  = ONE / (ONE + b*t2)

     h(i) = g3*gamma*log(ONE+delt*t2*q2)

     if (lpot) then

!    --------------------------------------------------------
!    energy done. now the potential, using appendix e of [b].
!    --------------------------------------------------------

       b2     = b*b
       g4     = g3*g91(i)
       t4     = t2*t2
       t6     = t4*t2
       rsthrd = rs(i)/THREE
       gz     = (((ONE+zet(i))**TWO+eta)**SIXTHM-((ONE-zet(i))**TWO+eta)**SIXTHM)/THREE
       fac    = delt/b + ONE
       bg     = -THREE*b2*ec(i)*fac/(bet*g4)
       bec    = b2*fac/(bet*g3)

!      ------------------------------------ NEW for MUcTWO ------------begin-
       q4     = ONE / (ONE + b*t2 + delt*t2)
       q2sq   = q2*q2
       q4sq   = q4*q4
       hb     = -gamma*g3*delt*t4*q2*q4
!      ------------------------------------ NEW for MUcTWO --------------end-

       hrs    = - rsthrd*hb*bec*ecrs(i)

!      ------------------------------------ NEW for MUcTWO ------------begin-
       fact0  = FOUR + FOUR*b*t2 + TWO*delt*t2
       hbt    = -gamma*g3*delt*t2*fact0*q2sq*q4sq
!      ------------------------------------ NEW for MUcTWO --------------end-

       hrst   = rsthrd*t2*hbt*bec*ecrs(i)
       hz     = THREE*gz*h(i)/g91(i) + hb*(bg*gz+bec*eczet(i))

!      ------------------------------------ NEW for MUcTWO ------------begin-
       ht     = gamma*g3*TWO*delt*q2*q4
!      ------------------------------------ NEW for MUcTWO --------------end-

       hzt    = THREE*gz*ht/g91(i) + hbt*(bg*gz+bec*eczet(i))
       fact2  = q4*q5 + b*t2*(q4*q9+q5)
       fact3  = TWO*b*q5*q9 + delt*fact2

!      ------------------------------------ NEW for MUcTWO ------------begin-
       fact1  = TWO*b+delt+TWO*t2*(b2+b*delt)
       htt    = -gamma*g3*FOUR*delt*t(i)*fact1*q2sq*q4sq
!      ------------------------------------ NEW for MUcTWO --------------end-

       comm   = h(i) + hrs + hrst + t2*ht/SIX + SEVEN*t2*t(i)*htt/SIX
       pref   = hz - gz*t2*ht/g91(i)
       fact5  = gz*(TWO*ht+t(i)*htt)/g91(i)
       comm   = comm - pref*zet(i) - uu(i)*htt - vv(i)*ht - ww(i)*(hzt-fact5)
       dvcup(i)  = comm + pref
       dvcdn(i)  = comm - pref

     else

       dvcup(i) = ZERO
       dvcdn(i) = ZERO

     endif

   enddo lpi

end subroutine corggapbe2