#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef HAVE_LIBXC
#include "xc_version.h"
#endif

!!@LICENSE
!
!******************************************************************************
! MODULE m_ldaxc
! Provides routines for LDA XC functional evaluation
!******************************************************************************
!
!   PUBLIC procedures available from this module:
! ldaxc   ! General subroutine for all coded LDA XC functionals
! exchng  ! Local exchange
! pw92c   ! Perdew & Wang, PRB, 45, 13244 (1992) (Correlation only)
! pw92xc  ! Perdew & Wang, PRB, 45, 13244 (1992)
! pzxc    ! Perdew & Zunger, PRB 23, 5075 (1981)
!
!   PUBLIC parameters, types, and variables available from this module:
! none
!
!******************************************************************************
!
!   USED module procedures:
! use gridxc_sys,        only: die             ! Termination routine
!
!   USED module parameters:
! use gridxc_precision,  only: dp              ! Real double precision type
!
!   EXTERNAL procedures used:
! none
!
!******************************************************************************

module gridxc_lda

  implicit none

  private

  public :: ldaxc   ! General subroutine for all coded LDA XC functionals
  public :: exchng  ! Local exchange
  public :: pw92c   ! Perdew & Wang, PRB, 45, 13244 (1992) (Correlation only)
  public :: pw92xc  ! Perdew & Wang, PRB, 45, 13244 (1992)
  public :: pzxc    ! Perdew & Zunger, PRB 23, 5075 (1981)

contains

  SUBROUTINE LDAXC( AUTHOR, IREL, nspin, D, EPSX, EPSC, VX, VC, &
      DVXDN, DVCDN &
#ifdef HAVE_LIBXC
      , is_libxc, xc_func, xc_info )
#else
    )
#endif

    ! ******************************************************************
    ! Finds the exchange and correlation energies and potentials, in the
    ! Local (spin) Density Approximation.
    ! Written by L.C.Balbas and J.M.Soler, Dec'96.
    ! Non-collinear spin added by J.M.Soler, May'98
    ! Interface with LibXC added by Micael Oliveira, Jul.2014
    ! *********** INPUT ************************************************
    ! CHARACTER*(*) AUTHOR : Parametrization desired:
    !     'CA' or 'PZ' => LSD Perdew & Zunger, PRB 23, 5075 (1981)
    !     'PW92' => LSD Perdew & Wang, PRB, 45, 13244 (1992)
    !     A check is made first for libxc identifiers (strings of the
    !     form "GGA_C_PBE", which are morphologically identical to the
    !     functional id macros but in string form.
    !     If the check fails the routine falls back on the built-in cases
    !     CA/PZ or PW92.
    !                     Uppercase is optional
    ! INTEGER IREL     : Relativistic exchange? (0=>no, 1=>yes)
    ! INTEGER nspin    : nspin=1 => unpolarized; nspin=2 => polarized;
    !                    nspin=4 => non-collinear polarization
    ! REAL*8  D(nspin) : Local (spin) density. For non-collinear
    !                    polarization, the density matrix is given by:
    !                    D(1)=D11, D(2)=D22, D(3)=Real(D12), D(4)=Im(D12)
    !
    !     Optional arguments if Libxc is compiled in (-DLIBXC):
    !     logical :: is_libxc    ! flag that determines whether it is a libxc functional
    !     type(xc_f90_pointer_t) :: xc_func, xc_info ! pre-initialized libxc objects
    
    !*********** OUTPUT ***********************************************
    ! REAL*8 EPSX, EPSC : Exchange and correlation energy densities
    ! REAL*8 VX(nspin), VC(nspin) : Exchange and correlation potentials,
    !                               defined as dExc/dD(ispin)
    ! REAL*8 DVXDN(nspin,nspin)  :  Derivative of exchange potential with
    !                               respect the charge density, defined 
    !                               as DVx(spin1)/Dn(spin2)
    ! REAL*8 DVCDN(nspin,nspin)  :  Derivative of correlation potential
    !                               respect the charge density, defined 
    !                               as DVc(spin1)/Dn(spin2)
    ! *********** UNITS ************************************************
    ! Lengths in Bohr, energies in Hartrees
    ! ******************************************************************
    
    use gridxc_precision, only : dp
    use gridxc_sys,       only : die

#ifdef HAVE_LIBXC
#if XC_MAJOR_VERSION >= 4
    use xc_f03_lib_m
    use iso_c_binding, only: c_size_t, c_int
#else
    use xc_f90_types_m
    use xc_f90_lib_m
#endif
#endif /* HAVE_LIBXC */

    CHARACTER*(*),intent(in) :: AUTHOR ! LDA flavour ('PZ'|'PW92')
    INTEGER, intent(in) :: IREL        ! Relativistic exchange? 0=>no, 1=>yes
    INTEGER, intent(in) :: nspin       ! Number of spin components
    real(dp),intent(in) :: D(nspin)    ! Electron density (matrix)
    real(dp),intent(out):: EPSX        ! Exchange energy per electron
    real(dp),intent(out):: EPSC        ! Correlation energy per electron
    real(dp),intent(out):: VX(nspin)   ! Exchange potential
    real(dp),intent(out):: VC(nspin)   ! Correlation potential
    real(dp),intent(out):: DVXDN(nspin,nspin) ! dVX(spin1)/dDens(spin2)
    real(dp),intent(out):: DVCDN(nspin,nspin) ! dVC(spin1)/dDens(spin2)
#ifdef HAVE_LIBXC
    logical, intent(in), optional    :: is_libxc
#if XC_MAJOR_VERSION >= 4
    type(xc_f03_func_t), optional :: xc_func
    type(xc_f03_func_info_t), optional :: xc_info
# if XC_MAJOR_VERSION >= 5
    integer(c_size_t), parameter :: one = 1
# else
    integer(c_int), parameter :: one = 1
# endif
#else
    type(xc_f90_pointer_t), optional :: xc_func, xc_info
#endif
#endif

    INTEGER           IS, NS, ISPIN1, ISPIN2
    real(dp)          DD(2), DPOL, DTOT, TINY, VCD(2), VPOL, VXD(2)

#ifdef HAVE_LIBXC
    logical  :: lis_libxc
    integer :: xc_kind
    real(dp) :: eps(1), ldedn(nspin), fxc(3), kxc(4), v4xc(5)
#endif

    PARAMETER ( TINY = 1.D-12 )

    IF (nspin .EQ. 4) THEN
      ! Find eigenvalues of density matrix (up and down densities
      ! along the spin direction)
      ! Note: D(1)=D11, D(2)=D22, D(3)=Real(D12), D(4)=Im(D12)
      NS = 2
      DTOT = D(1) + D(2)
      DPOL = SQRT( (D(1)-D(2))**2 + 4.D0*(D(3)**2+D(4)**2) )
      DD(1) = 0.5D0 * ( DTOT + DPOL )
      DD(2) = 0.5D0 * ( DTOT - DPOL )
    ELSE
      NS = nspin
      DO IS = 1,nspin
        !ag       Avoid negative densities
        DD(IS) = max(D(IS),0.0d0)
      end DO
    ENDIF


    DO ISPIN2 = 1, nspin
      DO ISPIN1 = 1, nspin
        DVXDN(ISPIN1,ISPIN2) = 0.D0
        DVCDN(ISPIN1,ISPIN2) = 0.D0
      ENDDO
    ENDDO

#ifdef HAVE_LIBXC
    if ( present(is_libxc) ) then
      lis_libxc = is_libxc
      if ( lis_libxc ) then
        if ((.not. present(xc_func)) .or. &
            (.not. present(xc_info))) then
          call die("xc_func and xc_info not present")
        endif
      endif
    else
      lis_libxc = .false.
    end if

    if (lis_libxc) then

      !     pass the arrays as assumed-size
      ! '1' is the number of grid points: the current approach is
      ! wasteful since it initializes and cleans the functional object
      ! at every grid point.
      !      call xc_f90_lda_exc_vxc(xc_func, 1, DD(1), eps, ldedn(1))
      ! Use general LDA libxc call to get fxc
#if XC_MAJOR_VERSION >= 4
      IF (xc_f03_func_info_get_family(xc_info) /= XC_FAMILY_LDA) THEN
        call die('LDAXC: Functional is not an LDA')
      ENDIF
      xc_kind = xc_f03_func_info_get_kind(xc_info)
# if XC_MAJOR_VERSION >= 5
      call xc_f03_lda_exc_vxc_fxc( xc_func, one, DD, eps, ldedn, fxc)
# else
      call xc_f03_lda( xc_func, one, DD, eps, ldedn, fxc, kxc)
# endif
#else
      IF (xc_f90_info_family(xc_info) /= XC_FAMILY_LDA) THEN
        call die('LDAXC: Functional is not an LDA')
      ENDIF
      xc_kind = xc_f90_info_kind(xc_info)
      call xc_f90_lda( xc_func, 1, DD(1), eps(1), ldedn(1), fxc(1), kxc(1))
#endif
      IF (xc_kind == XC_CORRELATION) THEN
        EPSC = eps(1)
        VCD(1:nspin) = ldedn(1:nspin)
        EPSX = 0.0_dp
        VXD(1:nspin) = 0.0_dp
        DVXDN(:,:) = 0.0_dp
        if (nspin == 1) then
           DVCDN(1,1) = fxc(1)
        else
           DVCDN(1,1) = fxc(1)
           DVCDN(2,2) = fxc(3)
           DVCDN(1,2) = fxc(2)
           DVCDN(2,1) = fxc(2)
        endif
      ELSE if (xc_kind == XC_EXCHANGE) THEN
        EPSX = eps(1)
        VXD(1:nspin) = ldedn(1:nspin)
        EPSC = 0.0_dp
        VCD(1:nspin) = 0.0_dp
        DVCDN(:,:) = 0.0_dp
        if (nspin == 1) then
           DVXDN(1,1) = fxc(1)
        else
           DVXDN(1,1) = fxc(1)
           DVXDN(2,2) = fxc(3)
           DVXDN(1,2) = fxc(2)
           DVXDN(2,1) = fxc(2)
        endif
      ELSE  ! combined functional, use an arbitrary 50/50 split
        EPSX = 0.5_dp * eps(1)
        EPSC = 0.5_dp * eps(1)
        VXD(1:nspin) = 0.5_dp * ldedn(1:nspin)
        VCD(1:nspin) = 0.5_dp * ldedn(1:nspin)
        if (nspin == 1) then
           DVXDN(1,1) = 0.5_dp * fxc(1)
        else
           DVXDN(1,1) = 0.5_dp * fxc(1)
           DVXDN(2,2) = 0.5_dp * fxc(3)
           DVXDN(1,2) = 0.5_dp * fxc(2)
           DVXDN(2,1) = 0.5_dp * fxc(2)
        endif
        DVCDN(:,:) = DVXDN(:,:)
      ENDIF

    ELSE IF ( AUTHOR.EQ.'CA' .OR. AUTHOR.EQ.'ca' .OR. &
#else
    IF ( AUTHOR.EQ.'CA' .OR. AUTHOR.EQ.'ca' .OR. &
#endif
      AUTHOR.EQ.'PZ' .OR. AUTHOR.EQ.'pz') THEN
      CALL PZXC( IREL, NS, DD, EPSX, EPSC, VXD, VCD, DVXDN, DVCDN )

    ELSEIF ( AUTHOR.EQ.'PW92' .OR. AUTHOR.EQ.'pw92' ) THEN
       ! Note that DV{XC}DN is not implemented for this flavor.
       ! Better use libxc
       CALL PW92XC( IREL, NS, DD, EPSX, EPSC, VXD, VCD )

    ELSE
      call die('LDAXC: Unknown author ' // trim(AUTHOR))
    ENDIF

    IF (nspin .EQ. 4) THEN
      ! NOTE: dE2/DN2 is not back-converted !
      ! Find dE/dD(ispin) = dE/dDup * dDup/dD(ispin) +
      !                     dE/dDdown * dDown/dD(ispin)
      VPOL  = (VXD(1)-VXD(2)) * (D(1)-D(2)) / (DPOL+TINY)
      VX(1) = 0.5D0 * ( VXD(1) + VXD(2) + VPOL )
      VX(2) = 0.5D0 * ( VXD(1) + VXD(2) - VPOL )
      VX(3) = (VXD(1)-VXD(2)) * D(3) / (DPOL+TINY)
      VX(4) = (VXD(1)-VXD(2)) * D(4) / (DPOL+TINY)
      VPOL  = (VCD(1)-VCD(2)) * (D(1)-D(2)) / (DPOL+TINY)
      VC(1) = 0.5D0 * ( VCD(1) + VCD(2) + VPOL )
      VC(2) = 0.5D0 * ( VCD(1) + VCD(2) - VPOL )
      VC(3) = (VCD(1)-VCD(2)) * D(3) / (DPOL+TINY)
      VC(4) = (VCD(1)-VCD(2)) * D(4) / (DPOL+TINY)
    ELSE
      DO IS = 1,nspin
        VX(IS) = VXD(IS)
        VC(IS) = VCD(IS)
      end do
    ENDIF
  END SUBROUTINE LDAXC


  subroutine exchng( IREL, NSP, DS, EX, VX )

    ! *****************************************************************
    !  Finds local exchange energy density and potential
    !  Adapted by J.M.Soler from routine velect of Froyen's 
    !    pseudopotential generation program. Madrid, Jan'97. Version 0.5.
    !  Relativistic exchange modified by JMS, May.2014
    ! **** Input ******************************************************
    ! INTEGER IREL    : relativistic-exchange switch (0=no, 1=yes)
    ! INTEGER NSP     : spin-polarizations (1=>unpolarized, 2=>polarized)
    ! REAL*8  DS(NSP) : total (nsp=1) or spin (nsp=2) electron density
    ! **** Output *****************************************************
    ! REAL*8  EX      : exchange energy density
    ! REAL*8  VX(NSP) : (spin-dependent) exchange potential
    ! **** Units ******************************************************
    ! Densities in electrons/Bohr**3
    ! Energies in Hartrees
    ! *****************************************************************

    use gridxc_precision, only: dp
    implicit none

    integer, intent(in) :: nsp, irel
    real(dp), intent(in)             :: DS(NSP)
    real(dp), intent(out)            :: VX(NSP)
    real(dp), intent(out)            :: EX

    real(dp), parameter :: zero = 0.0_dp, one = 1.0_dp
    real(dp), parameter :: pfive = 0.5_dp, opf = 1.5_dp
    !      real(dp), parameter :: c = 137.035999_dp       ! speed of light in a.u.
    !      real(dp), parameter :: C014 = 0.014_dp         ! (9*pi/4)^(1/3)/c
    real(dp), parameter :: C014 = 0.0140047747_dp  ! updated JMS, May.2014

    real(dp) :: a0, alp, sb, rs
    real(dp) :: pi, trd, ftrd, tftm
    real(dp) :: d1, d2, d, z, fz, fzp, vxp, exp_var
    real(dp) :: beta, vxf, exf, alb
    real(dp) :: dBETAdD, dETAdD, dGAMMAdD, dPHIdD, dKFdD
    real(dp) :: ETA, GAMMA, KF, PHI

    PI=4*ATAN(ONE)
    TRD = ONE/3
    FTRD = 4*TRD
    TFTM = 2._dp**FTRD-2

    IF (NSP .EQ. 2) THEN
      D1 = MAX(DS(1),ZERO)
      D2 = MAX(DS(2),ZERO)
      D = D1 + D2
      IF (D .LE. ZERO) THEN
        EX = ZERO
        VX(1) = ZERO
        VX(2) = ZERO
        RETURN
      ENDIF
      Z = (D1 - D2) / D
      FZ = ((1+Z)**FTRD+(1-Z)**FTRD-2)/TFTM
      FZP = FTRD*((1+Z)**TRD-(1-Z)**TRD)/TFTM 
    ELSE
      D = DS(1)
      IF (D .LE. ZERO) THEN
        EX = ZERO
        VX(1) = ZERO
        RETURN
      ENDIF
      Z = ZERO
      FZ = ZERO
      FZP = ZERO
    ENDIF

    A0 = (4/(9*PI))**TRD 
    ALP = 2 * TRD                      ! X-alpha parameter
    RS = (3 / (4*PI*D) )**TRD
    VXP = -(3*ALP/(2*PI*A0*RS))        ! VX=-KF/PI
    EXP_VAR = 3*VXP/4                  ! epsX=-3*KF/4/PI
    IF (IREL .EQ. 1) THEN
      ! Ref: MacDonald and Vosco, J.Phys C 12, 2977 (1979)
      BETA = C014/RS                   ! ratio of Fermi to light speed
      SB = SQRT(1+BETA*BETA)
      ALB = LOG(BETA+SB)
      VXP = VXP * (-PFIVE + OPF * ALB / (BETA*SB))
      EXP_VAR = EXP_VAR * (ONE-OPF*((BETA*SB-ALB)/BETA**2)**2) 
    ENDIF

    IF (NSP .EQ. 2) THEN
      VXF = 2**TRD*VXP
      EXF = 2**TRD*EXP_VAR
      VX(1) = VXP + FZ*(VXF-VXP) + (1-Z)*FZP*(EXF-EXP_VAR)
      VX(2) = VXP + FZ*(VXF-VXP) - (1+Z)*FZP*(EXF-EXP_VAR)
      EX    = EXP_VAR + FZ*(EXF-EXP_VAR)
    ELSE
      VX(1) = VXP
      EX    = EXP_VAR
    ENDIF
  END subroutine exchng



  SUBROUTINE PW92C( nspin, Dens, EC, VC )

    ! ********************************************************************
    ! Implements the Perdew-Wang'92 local correlation (beyond RPA).
    ! Ref: J.P.Perdew & Y.Wang, PRB, 45, 13244 (1992)
    ! Written by L.C.Balbas and J.M.Soler. Dec'96.  Version 0.5.
    ! ********* INPUT ****************************************************
    ! INTEGER nspin       : Number of spin polarizations (1 or 2)
    ! REAL*8  Dens(nspin) : Local (spin) density
    ! ********* OUTPUT ***************************************************
    ! REAL*8  EC        : Correlation energy density
    ! REAL*8  VC(nspin) : Correlation (spin) potential
    ! ********* UNITS ****************************************************
    ! Densities in electrons per Bohr**3
    ! Energies in Hartrees
    ! ********* ROUTINES CALLED ******************************************
    ! None
    ! ********************************************************************

    use gridxc_precision, only : dp

    ! Next line is nonstandard but may be supressed
    implicit          none

    ! Argument types and dimensions
    INTEGER           nspin
    real(dp)          Dens(nspin), EC, VC(nspin)

    ! Internal variable declarations
    INTEGER           IG
    real(dp)          A(0:2), ALPHA1(0:2), B, BETA(0:2,4), C, &
                      DBDRS, DECDD(2), DECDRS, DECDZ, DENMIN, DFDZ, &
                      DGDRS(0:2), DCDRS, DRSDD, DTOT, DZDD(2), &
                      F, FPP0, FOUTHD, G(0:2), HALF, ONE, &
                     P(0:2), PI, RS, THD, THRHLF, ZETA

    ! Add tiny numbers to avoid numerical errors
    PARAMETER ( DENMIN = 1.D-12 )
    PARAMETER ( ONE    = 1.D0 + 1.D-12 )

    ! Fix some numerical constants
    PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, &
                    THD=1.D0/3.D0, THRHLF=1.5D0 )

    ! Parameters from Table I of Perdew & Wang, PRB, 45, 13244 (92)
    DATA P      / 1.00d0,     1.00d0,     1.00d0     /
    DATA A      / 0.031091d0, 0.015545d0, 0.016887d0 /
    DATA ALPHA1 / 0.21370d0,  0.20548d0,  0.11125d0  /
    DATA BETA   / 7.5957d0,  14.1189d0,  10.357d0, &
                      3.5876d0,   6.1977d0,   3.6231d0, &
                  1.6382d0,   3.3662d0,   0.88026d0, &
                  0.49294d0,  0.62517d0,  0.49671d0 /

    ! Find rs and zeta
    PI = 4 * ATAN(1.D0)
    IF (nspin .EQ. 1) THEN
      DTOT = MAX( DENMIN, Dens(1) )
      ZETA = 0
      RS = ( 3 / (4*PI*DTOT) )**THD
      !       Find derivatives dRs/dDens and dZeta/dDens
      DRSDD = (- RS) / DTOT / 3
      DZDD(1) = 0
    ELSE
      DTOT = MAX( DENMIN, Dens(1)+Dens(2) )
      ZETA = ( Dens(1) - Dens(2) ) / DTOT
      RS = ( 3 / (4*PI*DTOT) )**THD
      DRSDD = (- RS) / DTOT / 3
      DZDD(1) =   (ONE - ZETA) / DTOT
      DZDD(2) = - (ONE + ZETA) / DTOT
    ENDIF

    ! Find eps_c(rs,0)=G(0), eps_c(rs,1)=G(1) and -alpha_c(rs)=G(2)
    ! using eq.(10) of cited reference (Perdew & Wang, PRB, 45, 13244 (92))
    DO IG = 0,2
      B = BETA(IG,1) * RS**HALF   + &
          BETA(IG,2) * RS         + &
          BETA(IG,3) * RS**THRHLF + &
          BETA(IG,4) * RS**(P(IG)+1)
      DBDRS = BETA(IG,1) * HALF      / RS**HALF + &
          BETA(IG,2)                         + &
          BETA(IG,3) * THRHLF    * RS**HALF + &
          BETA(IG,4) * (P(IG)+1) * RS**P(IG)
      C = 1 + 1 / (2 * A(IG) * B)
      DCDRS = - ( (C-1) * DBDRS / B )
      G(IG) = (- 2) * A(IG) * ( 1 + ALPHA1(IG)*RS ) * LOG(C)
      DGDRS(IG) = (- 2) *A(IG) * ( ALPHA1(IG) * LOG(C) + &
          (1+ALPHA1(IG)*RS) * DCDRS / C )
    end do

    ! Find f''(0) and f(zeta) from eq.(9)
    C = 1 / (2**FOUTHD - 2)
    FPP0 = 8 * C / 9
    F = ( (ONE+ZETA)**FOUTHD + (ONE-ZETA)**FOUTHD - 2 ) * C
    DFDZ = FOUTHD * ( (ONE+ZETA)**THD - (ONE-ZETA)**THD ) * C

    ! Find eps_c(rs,zeta) from eq.(8)
    EC = G(0) - G(2) * F / FPP0 * (ONE-ZETA**4) + &
        (G(1)-G(0)) * F * ZETA**4
    DECDRS = DGDRS(0) - DGDRS(2) * F / FPP0 * (ONE-ZETA**4) + &
        (DGDRS(1)-DGDRS(0)) * F * ZETA**4
    DECDZ = (- G(2)) / FPP0 * ( DFDZ*(ONE-ZETA**4) - F*4*ZETA**3 ) + &
        (G(1)-G(0)) * ( DFDZ*ZETA**4 + F*4*ZETA**3 )

    ! Find correlation potential
    IF (nspin .EQ. 1) THEN
      DECDD(1) = DECDRS * DRSDD
      VC(1) = EC + DTOT * DECDD(1)
    ELSE
      DECDD(1) = DECDRS * DRSDD + DECDZ * DZDD(1)
      DECDD(2) = DECDRS * DRSDD + DECDZ * DZDD(2)
      VC(1) = EC + DTOT * DECDD(1)
      VC(2) = EC + DTOT * DECDD(2)
    ENDIF

  END SUBROUTINE PW92C



  SUBROUTINE PW92XC( IREL, nspin, Dens, EPSX, EPSC, VX, VC )

    ! ********************************************************************
    ! Implements the Perdew-Wang'92 LDA/LSD exchange correlation
    ! Ref: J.P.Perdew & Y.Wang, PRB, 45, 13244 (1992)
    ! Written by L.C.Balbas and J.M.Soler. Dec'96. Version 0.5.
    ! ********* INPUT ****************************************************
    ! INTEGER IREL        : Relativistic-exchange switch (0=No, 1=Yes)
    ! INTEGER nspin       : Number of spin polarizations (1 or 2)
    ! REAL*8  Dens(nspin) : Local (spin) density
    ! ********* OUTPUT ***************************************************
    ! REAL*8  EPSX       : Exchange energy density
    ! REAL*8  EPSC       : Correlation energy density
    ! REAL*8  VX(nspin)  : Exchange (spin) potential
    ! REAL*8  VC(nspin)  : Correlation (spin) potential
    ! ********* UNITS ****************************************************
    ! Densities in electrons per Bohr**3
    ! Energies in Hartrees
    ! ********* ROUTINES CALLED ******************************************
    ! EXCHNG, PW92C
    ! ********************************************************************

    use gridxc_precision, only : dp

    INTEGER           IREL, nspin
    real(dp)          Dens(nspin), EPSX, EPSC, VC(nspin), VX(nspin)

    CALL EXCHNG( IREL, nspin, Dens, EPSX, VX )
    CALL PW92C( nspin, Dens, EPSC, VC )
    
  END SUBROUTINE PW92XC



  SUBROUTINE PZXC( IREL, NSP, DS, EX, EC, VX, VC, DVXDN, DVCDN )

    ! *****************************************************************
    !  Perdew-Zunger parameterization of Ceperley-Alder exchange and 
    !  correlation. Ref: Perdew & Zunger, Phys. Rev. B 23 5075 (1981).
    !  Adapted by J.M.Soler from routine velect of Froyen's 
    !    pseudopotential generation program. Madrid, Jan'97.
    ! **** Input *****************************************************
    ! INTEGER IREL    : relativistic-exchange switch (0=no, 1=yes)
    ! INTEGER NSP     : spin-polarizations (1=>unpolarized, 2=>polarized)
    ! REAL*8  DS(NSP) : total (nsp=1) or spin (nsp=2) electron density
    ! **** Output *****************************************************
    ! REAL*8  EX            : exchange energy density
    ! REAL*8  EC            : correlation energy density
    ! REAL*8  VX(NSP)       : (spin-dependent) exchange potential
    ! REAL*8  VC(NSP)       : (spin-dependent) correlation potential
    ! REAL*8  DVXDN(NSP,NSP): Derivative of the exchange potential
    !                         respect the charge density, 
    !                         Dvx(spin1)/Dn(spin2)
    ! REAL*8  DVCDN(NSP,NSP): Derivative of the correlation potential
    !                         respect the charge density, 
    !                         Dvc(spin1)/Dn(spin2)
    ! **** Units *******************************************************
    ! Densities in electrons/Bohr**3
    ! Energies in Hartrees
    ! *****************************************************************

    use gridxc_precision, only: dp

    implicit none

    integer  :: nsp, irel, isp1, isp2, isp
    real(dp) :: DS(NSP), VX(NSP), VC(NSP), &
        DVXDN(NSP,NSP), DVCDN(NSP,NSP)
    real(dp), parameter :: &
        ZERO=0.D0,ONE=1.D0,PFIVE=.5D0,OPF=1.5D0,PNN=.99D0, &
        PTHREE=0.3D0,PSEVF=0.75D0,C0504=0.0504D0, &
        C0254=0.0254D0,C014=0.014D0,C0406=0.0406D0, &
        C15P9=15.9D0,C0666=0.0666D0,C11P4=11.4D0, &
        C045=0.045D0,C7P8=7.8D0,C88=0.88D0,C20P59=20.592D0, &
        C3P52=3.52D0,C0311=0.0311D0,C0014=0.0014D0, &
        C0538=0.0538D0,C0096=0.0096D0,C096=0.096D0, &
        C0622=0.0622D0,C004=0.004D0,C0232=0.0232D0, &
        C1686=0.1686D0,C1P398=1.3981D0,C2611=0.2611D0, &
        C2846=0.2846D0,C1P053=1.0529D0,C3334=0.3334D0

    !    Ceperly-Alder 'ca' constants. Internal energies in Rydbergs.
    real(dp), parameter :: &
        CON1=1.D0/6, CON2=0.008D0/3, CON3=0.3502D0/3, &
        CON4=0.0504D0/3, CON5=0.0028D0/3, CON6=0.1925D0/3, &
        CON7=0.0206D0/3, CON8=9.7867D0/6, CON9=1.0444D0/3, &
        CON10=7.3703D0/6, CON11=1.3336D0/3

    !      X-alpha parameter:
    real(dp), PARAMETER :: ALP = 2.D0 / 3.D0 

    !      Other variables converted into parameters by J.M.Soler
    real(dp), parameter :: &
        TINY = 1.D-6 , &
        PI   = 3.14159265358979312_dp, &
        TWO  = 2.0D0, &
        HALF = 0.5D0, &
        TRD  = 1.D0 / 3.D0, &
        FTRD = 4.D0 / 3.D0, &
        TFTM = 0.51984209978974638D0, &
        A0   = 0.52106176119784808D0, &
        CRS  = 0.620350490899400087D0, &
        CXP  = (- 3.D0) * ALP / (PI*A0), &
        CXF  = 1.25992104989487319D0 

    real(dp)  :: d1, d2, d, z, fz, fzp
    real(dp)  :: ex, ec, dfzpdn, rs, vxp, exp_var
    real(dp)  :: beta, sb, alb, vxf, exf, dvxpdn
    real(dp)  :: dvxfdn, sqrs, te, be, ecp, vcp
    real(dp)  :: dtedn, be2, dbedn, dvcpdn, decpdn
    real(dp)  :: ecf, vcf, dvcfdn, decfdn, rslog


    !      Find density and polarization
    IF (NSP .EQ. 2) THEN
      D1 = MAX(DS(1),ZERO)
      D2 = MAX(DS(2),ZERO)
      D = D1 + D2
      IF (D .LE. ZERO) THEN
        EX = ZERO
        EC = ZERO
        VX(1) = ZERO
        VX(2) = ZERO
        VC(1) = ZERO
        VC(2) = ZERO
        RETURN
      ENDIF
      !
      !        Robustness enhancement by Jose Soler (August 2002)
      !
      Z = (D1 - D2) / D
      IF (Z .LE. -ONE) THEN
        FZ = (TWO**FTRD-TWO)/TFTM
        FZP = -FTRD*TWO**TRD/TFTM
        DFZPDN = FTRD*TRD*TWO**(-ALP)/TFTM
      ELSEIF (Z .GE. ONE) THEN
        FZ = (TWO**FTRD-TWO)/TFTM
        FZP = FTRD*TWO**TRD/TFTM
        DFZPDN = FTRD*TRD*TWO**(-ALP)/TFTM
      ELSE
        FZ = ((ONE+Z)**FTRD+(ONE-Z)**FTRD-TWO)/TFTM
        FZP = FTRD*((ONE+Z)**TRD-(ONE-Z)**TRD)/TFTM 
        DFZPDN = FTRD*TRD*((ONE+Z)**(-ALP) + (ONE-Z)**(-ALP))/TFTM
      ENDIF
    ELSE
      D = DS(1)
      IF (D .LE. ZERO) THEN
        EX = ZERO
        EC = ZERO
        VX(1) = ZERO
        VC(1) = ZERO
        RETURN
      ENDIF
      Z = ZERO
      FZ = ZERO
      FZP = ZERO
    ENDIF
    RS = CRS / D**TRD

    !      Exchange
    VXP = CXP / RS
    EXP_VAR = 0.75D0 * VXP
    IF (IREL .EQ. 1) THEN
      BETA = C014/RS
      IF (BETA .LT. TINY) THEN
        SB = ONE + HALF*BETA**2
        ALB = BETA
      ELSE
        SB = SQRT(1+BETA*BETA)
        ALB = LOG(BETA+SB)
      ENDIF
      VXP = VXP * (-PFIVE + OPF * ALB / (BETA*SB))
      EXP_VAR = EXP_VAR *(ONE-OPF*((BETA*SB-ALB)/BETA**2)**2) 
    ENDIF
    VXF = CXF * VXP
    EXF = CXF * EXP_VAR
    DVXPDN = TRD * VXP / D
    DVXFDN = TRD * VXF / D

    !      Correlation 
    IF (RS .GT. ONE) THEN  
      SQRS=SQRT(RS)
      TE = ONE+CON10*SQRS+CON11*RS
      BE = ONE+C1P053*SQRS+C3334*RS
      ECP = -(C2846/BE)
      VCP = ECP*TE/BE
      DTEDN = ((CON10 * SQRS *HALF) + CON11 * RS)*(-TRD/D)
      BE2 = BE * BE
      DBEDN = ((C1P053 * SQRS *HALF) + C3334 * RS)*(-TRD/D)
      DVCPDN = -(C2846/BE2)*(DTEDN - 2.0D0 * TE * DBEDN/BE)
      DECPDN = (C2846/BE2)*DBEDN
      TE = ONE+CON8*SQRS+CON9*RS
      BE = ONE+C1P398*SQRS+C2611*RS
      ECF = -(C1686/BE)
      VCF = ECF*TE/BE
      DTEDN = ((CON8 * SQRS * HALF) + CON9 * RS)*(-TRD/D)
      BE2 = BE * BE
      DBEDN = ((C1P398 * SQRS * HALF) + C2611 * RS)*(-TRD/D)
      DVCFDN = -(C1686/BE2)*(DTEDN - 2.0D0 * TE * DBEDN/BE)
      DECFDN = (C1686/BE2)*DBEDN
    ELSE
      RSLOG=LOG(RS)
      ECP=(C0622+C004*RS)*RSLOG-C096-C0232*RS
      VCP=(C0622+CON2*RS)*RSLOG-CON3-CON4*RS
      DVCPDN = (CON2*RS*RSLOG + (CON2-CON4)*RS + C0622)*(-TRD/D)
      DECPDN = (C004*RS*RSLOG + (C004-C0232)*RS + C0622)*(-TRD/D)
      ECF=(C0311+C0014*RS)*RSLOG-C0538-C0096*RS
      VCF=(C0311+CON5*RS)*RSLOG-CON6-CON7*RS
      DVCFDN = (CON5*RS*RSLOG + (CON5-CON7)*RS + C0311)*(-TRD/D)
      DECFDN = (C0014*RS*RSLOG + (C0014-C0096)*RS + C0311)*(-TRD/D)
    ENDIF

    ISP1 = 1
    ISP2 = 2

    !      Find up and down potentials
    IF (NSP .EQ. 2) THEN
      EX    = EXP_VAR + FZ*(EXF-EXP_VAR)
      EC    = ECP + FZ*(ECF-ECP)
      VX(1) = VXP + FZ*(VXF-VXP) + (ONE-Z)*FZP*(EXF-EXP_VAR)
      VX(2) = VXP + FZ*(VXF-VXP) - (ONE+Z)*FZP*(EXF-EXP_VAR)
      VC(1) = VCP + FZ*(VCF-VCP) + (ONE-Z)*FZP*(ECF-ECP)
      VC(2) = VCP + FZ*(VCF-VCP) - (ONE+Z)*FZP*(ECF-ECP)

      !        Derivatives of exchange potential respect the density

      DVXDN(ISP1,ISP1) = &
          DVXPDN &
          +  FZP*(VXF-VXP-EXF+EXP_VAR)*( 2.D0*D2/(D*D) ) &
          +  FZ*(DVXFDN-DVXPDN)+(1-Z)*FZP*(VXF-VXP)/(4.D0*D) &
          +  (1-Z)*DFZPDN*(EXF-EXP_VAR)*( 2.D0*D2/(D*D) )
      DVXDN(ISP1,ISP2) = &
          DVXPDN &
          +  FZP*(VXF-VXP-EXF+EXP_VAR)*(-2.D0*D1/(D*D) ) &
          +  FZ*(DVXFDN-DVXPDN)+(1-Z)*FZP*(VXF-VXP)/(4.D0*D) &
          +  (1-Z)*DFZPDN*(EXF-EXP_VAR)*( -2.D0*D1/(D*D) )
      DVXDN(ISP2,ISP1) = &
          DVXPDN &
          +  FZP*(VXF-VXP-EXF+EXP_VAR)*( 2.D0*D2/(D*D) ) &
          +  FZ*(DVXFDN-DVXPDN)-(1+Z)*FZP*(VXF-VXP)/(4.D0*D) &
          -  (1+Z)*DFZPDN*(EXF-EXP_VAR)*( 2.D0*D2/(D*D) )
      DVXDN(ISP2,ISP2) = &
          DVXPDN &
          +  FZP*(VXF-VXP-EXF+EXP_VAR)*(-2.D0*D1/(D*D) ) &
          +  FZ*(DVXFDN-DVXPDN)-(1+Z)*FZP*(VXF-VXP)/(4.D0*D) &
          -  (1+Z)*DFZPDN*(EXF-EXP_VAR)*( -2.D0*D1/(D*D) )

      !        Derivatives of correlation potential respect the density

      DVCDN(ISP1,ISP1) = &
          DVCPDN &
          + FZP*(VCF-VCP-ECF+ECP)*( 2.D0*D2/(D*D) ) &
          + FZ*(DVCFDN-DVCPDN)+ (1-Z)*FZP*(DECFDN-DECPDN) &
          + (1-Z)*DFZPDN*(ECF-ECP)*( 2.D0*D2/(D*D) )
      DVCDN(ISP1,ISP2) = &
          DVCPDN &
          + FZP*(VCF-VCP-ECF+ECP)*(-2.D0*D1/(D*D) ) &
          + FZ*(DVCFDN-DVCPDN)+ (1-Z)*FZP*(DECFDN-DECPDN) &
          + (1-Z)*DFZPDN*(ECF-ECP)*( -2.D0*D1/(D*D) )
      DVCDN(ISP2,ISP1) = &
          DVCPDN &
          + FZP*(VCF-VCP-ECF+ECP)*( 2.D0*D2/(D*D) ) &
          + FZ*(DVCFDN-DVCPDN)- (1+Z)*FZP*(DECFDN-DECPDN) &
          - (1+Z)*DFZPDN*(ECF-ECP)*( 2.D0*D2/(D*D) )
      DVCDN(ISP2,ISP2) = &
          DVCPDN &
          + FZP*(VCF-VCP-ECF+ECP)*(-2.D0*D1/(D*D) ) &
          + FZ*(DVCFDN-DVCPDN)- (1+Z)*FZP*(DECFDN-DECPDN) &
          - (1+Z)*DFZPDN*(ECF-ECP)*( -2.D0*D1/(D*D) )

    ELSE
      EX    = EXP_VAR
      EC    = ECP
      VX(1) = VXP
      VC(1) = VCP
      DVXDN(1,1) = DVXPDN
      DVCDN(1,1) = DVCPDN
    ENDIF

    !      Change from Rydbergs to Hartrees
    EX = HALF * EX
    EC = HALF * EC
    DO ISP = 1,NSP
      VX(ISP) = HALF * VX(ISP)
      VC(ISP) = HALF * VC(ISP)
      DO ISP2 = 1,NSP
        DVXDN(ISP,ISP2) = HALF * DVXDN(ISP,ISP2)
        DVCDN(ISP,ISP2) = HALF * DVCDN(ISP,ISP2)
      end DO
    end DO
  END SUBROUTINE PZXC

END MODULE gridxc_lda

