Merge lp:~nickpapior/libgridxc/fortran-standard into lp:libgridxc

[Nick Papior]

I have now moved all codes to the fortran 90 standard.
There are, however still some goto statements which requires some major restructuring of the code. This could potentially be problematic.

Secondly, the m_fft_gpfa code is _very_ difficult to translate to f90 standard. Thus I have left it as is.

@Yann, I have requested your review since you know the limitations of build-systems.

@Alberto/@Jose, feel free to comment on this!

[Nick Papior (nickpapior)]

Note I have checked the output of test1 and test3 @ rev39 and my branch, they both yield _exactly_ the same output.

test3 has also been runned with mpirun -np 4.

[Yann Pouillon (pouillon)]

Great work! I know how tedious this kind of things are.

review: Approve

=== added directory 'docs/deprecated_routines'
=== added file 'docs/deprecated_routines/test_bessph.F90'
--- docs/deprecated_routines/test_bessph.F90	1970-01-01 00:00:00 +0000
+++ docs/deprecated_routines/test_bessph.F90	2018-02-19 19:59:41 +0000
@@ -0,0 +1,26 @@
+program test_bessph_update
+
+  ! One should get the bessph.f file from r39 and older
+  use m_bessph_old, only: bessph_old => bessph
+  ! This is the new bessph module
+  use m_bessph
+  
+  integer, parameter :: dp = kind(1.d0)
+  integer :: l, i
+  real(dp) :: x, dx, b1, b2
+  
+  dx = 0.00001_dp
+  do l = 0, 10
+    do i = 1, 10000000
+      x = (i-1) * dx
+      b1 = bessph_old(l, x)
+      b2 = bessph(l, x)
+      
+      if ( abs(b1 - b2) > 1.e-17_dp ) then
+        print *, 'Different: ', l, x, b1, b2
+      end if
+      
+    end do
+  end do
+  
+end program test_bessph_update

=== renamed file 'src/am05.f90' => 'src/am05.F90'
=== added file 'src/bessph.F90'
--- src/bessph.F90	1970-01-01 00:00:00 +0000
+++ src/bessph.F90	2018-02-19 19:59:41 +0000
@@ -0,0 +1,83 @@
+!!@LICENSE
+
+module m_bessph
+
+  implicit none
+
+  private
+
+  public :: bessph
+
+contains
+
+  function bessph(L,X) result(B)
+
+    ! RETURNS THE SPHERICAL BESSEL FUNCTION JL(X).
+    ! REF: ABRAMOWITZ AND STEGUN, FORMULAS 10.1.2 AND 10.1.19
+    ! WRITTEN BY J.SOLER. NOV/89.
+    ! F90 conversion, by Nick Papior, feb/2018
+    use precision, only: dp
+    use sys, only: die
+
+    integer, intent(in)  :: L
+    real(dp), intent(in) :: X
+    real(dp) :: B
+    
+    integer, parameter :: nterms = 100
+    real(dp), parameter :: zero = 0.0_dp, one = 1.0_dp, tiny=1.0e-15_dp
+
+    integer :: i, n
+    real(dp) :: switch, term, x2, y, fnm1, fn, fnp1
+    character(len=128):: msg
+
+    if ( abs(x) < max(1, 2*L - 1) ) then
+      ! Use power series
+      TERM = ONE
+      do I = 1 , L
+        TERM = TERM * X/(2*I+1)
+      end do
+      X2 = X * X
+      FN = ZERO
+      do I = 1 , NTERMS
+        FN = FN + TERM
+        TERM = -TERM * X2 / (4._dp * I * (I + L + 0.5_dp))
+        if ( abs(TERM) < tiny ) exit
+      end do
+      if ( abs(TERM) > tiny ) then
+        write(msg,*) 'BESSPH: SERIES HAS NOT CONVERGED. L,X=',L,X
+        call die(trim(msg))
+      end if
+      B = FN
+      
+    else
+      
+      ! Use explicit expressions or recurrence relation
+      select case ( L )
+      case ( 0 )
+        
+        B = sin(X) / X
+        
+      case ( 1 )
+        
+        B = (sin(X) / X - cos(X) ) / X
+        
+      case default
+        
+        Y = ONE / X
+        fnm1 = sin(X) * Y
+        fn = (fnm1 - cos(X)) * Y
+        DO N = 1 , L-1
+          FNP1 = (2 * N + 1) * Y * fn - fnm1
+          FNM1 = FN
+          FN = FNP1
+        end do
+        B = fn
+        
+      end select
+      
+    end if
+            
+  end function bessph
+  
+end module m_bessph
+

=== removed file 'src/bessph.f'
--- src/bessph.f	2016-05-09 12:15:45 +0000
+++ src/bessph.f	1970-01-01 00:00:00 +0000
@@ -1,66 +0,0 @@
-!!@LICENSE
-
-      MODULE m_bessph
-
-      CONTAINS
-
-      FUNCTION BESSPH (L,X)
-
-!  RETURNS THE SPHERICAL BESSEL FUNCTION JL(X).
-!  REF: ABRAMOWITZ AND STEGUN, FORMULAS 10.1.2 AND 10.1.19
-!  WRITTEN BY J.SOLER. NOV/89.
-
-      use precision, only: dp
-      use sys, only: die
-
-      integer, intent(in)  :: L
-      real(dp), intent(in) :: X
-      real(dp)             :: BESSPH
-
-      integer, parameter :: nterms = 100
-      real(dp),parameter ::
-     .   zero = 0.0_dp, one = 1.0_dp, tiny=1.0e-15_dp
-
-      integer :: i, n
-      real(dp) :: switch, term, x2, sum, y, fnm1, fn, fnp1
-      character(len=132):: msg
-
-      SWITCH=MAX(1,2*L-1)
-      IF (ABS(X).LT.SWITCH) THEN
-!        USE POWER SERIES
-         TERM=ONE
-         DO 10 I=1,L
-            TERM=TERM*X/(2*I+1)
-   10    CONTINUE
-         X2=X*X
-         SUM=ZERO
-         DO 20 I=1,NTERMS
-            SUM=SUM+TERM
-            TERM=(-TERM)*X2/(2*I*(2*I+2*L+1))
-            IF (ABS(TERM).LT.TINY) GO TO 30
-   20    CONTINUE
-            WRITE(msg,*) 'BESSPH: SERIES HAS NOT CONVERGED. L,X=',L,X
-            call die(trim(msg))
-   30    BESSPH=SUM
-      ELSE
-!        USE EXPLICIT EXPRESSIONS OR RECURRENCE RELATION
-         IF (L.EQ.0) THEN
-            BESSPH=SIN(X)/X
-         ELSEIF (L.EQ.1) THEN
-            BESSPH=(SIN(X)/X-COS(X))/X
-         ELSE
-            Y=ONE/X
-            FNM1=SIN(X)*Y
-            FN=(FNM1-COS(X))*Y
-            DO 40 N=1,L-1
-               FNP1=(2*N+1)*Y*FN-FNM1
-               FNM1=FN
-               FN=FNP1
-   40       CONTINUE
-            BESSPH=FN
-         ENDIF
-      ENDIF
-      END FUNCTION BESSPH
-
-      END MODULE m_bessph
-

=== renamed file 'src/cellsubs.f' => 'src/cellsubs.F90'
--- src/cellsubs.f	2016-05-09 12:15:45 +0000
+++ src/cellsubs.F90	2018-02-19 19:59:41 +0000
@@ -1,58 +1,64 @@
 !!@LICENSE
 
-      MODULE cellSubs
-
-      USE precision, only: dp   ! Double precision real kind
-
-      PUBLIC::  
-     .  reclat, ! Reciprocal unit cell vectors
-     .  volcel  ! Volume of the unit cell
-
-      PRIVATE
-
-      CONTAINS
-
-!--------------------------------------------------------------------
-
-      SUBROUTINE RECLAT (A,B,IOPT)
-
-C  CALCULATES RECIPROCAL LATTICE VECTORS. THEIR PRODUCT WITH DIRECT
-C  LATTICE VECTORS IS 1 IF IOPT=0 OR 2*PI IF IOPT=1
-
-      integer :: iopt, i
-      real(dp):: A(3,3),B(3,3), pi, c , ci
-      PI=ACOS(-1.D0)
-      B(1,1)=A(2,2)*A(3,3)-A(3,2)*A(2,3)
-      B(2,1)=A(3,2)*A(1,3)-A(1,2)*A(3,3)
-      B(3,1)=A(1,2)*A(2,3)-A(2,2)*A(1,3)
-      B(1,2)=A(2,3)*A(3,1)-A(3,3)*A(2,1)
-      B(2,2)=A(3,3)*A(1,1)-A(1,3)*A(3,1)
-      B(3,2)=A(1,3)*A(2,1)-A(2,3)*A(1,1)
-      B(1,3)=A(2,1)*A(3,2)-A(3,1)*A(2,2)
-      B(2,3)=A(3,1)*A(1,2)-A(1,1)*A(3,2)
-      B(3,3)=A(1,1)*A(2,2)-A(2,1)*A(1,2)
-      C=1.D0
-      IF (IOPT.EQ.1) C=2.D0*PI
-      DO 20 I=1,3
-         CI=C/(A(1,I)*B(1,I)+A(2,I)*B(2,I)+A(3,I)*B(3,I))
-         B(1,I)=B(1,I)*CI
-         B(2,I)=B(2,I)*CI
-         B(3,I)=B(3,I)*CI
-  20  CONTINUE
-      END SUBROUTINE RECLAT
-
-!--------------------------------------------------------------------
-
-      DOUBLE PRECISION FUNCTION VOLCEL( C )
-
-C  CALCULATES THE VOLUME OF THE UNIT CELL
-
-      real(dp) C(3,3)
-      VOLCEL = ( C(2,1)*C(3,2) - C(3,1)*C(2,2) ) * C(1,3) +
-     .         ( C(3,1)*C(1,2) - C(1,1)*C(3,2) ) * C(2,3) +
-     .         ( C(1,1)*C(2,2) - C(2,1)*C(1,2) ) * C(3,3)
-      VOLCEL = ABS( VOLCEL )
-      END FUNCTION VOLCEL
-
-      END MODULE cellSubs
+module cellSubs
+
+  use precision, only: dp   ! Double precision real kind
+
+  implicit none
+
+  private
+
+  public :: reclat ! Reciprocal unit cell vectors
+  public :: volcel  ! Volume of the unit cell
+
+contains
+
+  subroutine reclat(A,B,IOPT)
+    
+    !  CALCULATES RECIPROCAL LATTICE VECTORS. THEIR PRODUCT WITH DIRECT
+    !  LATTICE VECTORS IS 1 IF IOPT=0 OR 2*PI IF IOPT=1
+    real(dp), intent(in) :: A(3,3)
+    real(dp), intent(out) :: B(3,3)
+    integer, intent(in) :: iopt
+
+    real(dp), parameter :: PI = 3.14159265358979323846264338327950288419716939937510_dp
+
+    integer :: i
+    real(dp):: c , ci
+    
+    B(1,1) = A(2,2)*A(3,3)-A(3,2)*A(2,3)
+    B(2,1) = A(3,2)*A(1,3)-A(1,2)*A(3,3)
+    B(3,1) = A(1,2)*A(2,3)-A(2,2)*A(1,3)
+    B(1,2) = A(2,3)*A(3,1)-A(3,3)*A(2,1)
+    B(2,2) = A(3,3)*A(1,1)-A(1,3)*A(3,1)
+    B(3,2) = A(1,3)*A(2,1)-A(2,3)*A(1,1)
+    B(1,3) = A(2,1)*A(3,2)-A(3,1)*A(2,2)
+    B(2,3) = A(3,1)*A(1,2)-A(1,1)*A(3,2)
+    B(3,3) = A(1,1)*A(2,2)-A(2,1)*A(1,2)
+    C=1.D0
+    if ( IOPT == 1 ) C = 2.D0*PI
+    do I = 1 , 3
+      CI = C/(A(1,I)*B(1,I)+A(2,I)*B(2,I)+A(3,I)*B(3,I))
+      B(1,I) = B(1,I)*CI
+      B(2,I) = B(2,I)*CI
+      B(3,I) = B(3,I)*CI
+    end do
+    
+  end subroutine reclat
+
+  !--------------------------------------------------------------------
+
+  function volcel( C ) result(V)
+    !  CALCULATES THE VOLUME OF THE UNIT CELL
+    real(dp), intent(in) :: C(3,3)
+    real(dp) :: V
+    
+    V = ( C(2,1)*C(3,2) - C(3,1)*C(2,2) ) * C(1,3) + &
+        ( C(3,1)*C(1,2) - C(1,1)*C(3,2) ) * C(2,3) + &
+        ( C(1,1)*C(2,2) - C(2,1)*C(1,2) ) * C(3,3)
+    V = ABS( V )
+    
+  end function volcel
+
+end module cellsubs
 

=== renamed file 'src/chkgmx.f' => 'src/chkgmx.F90'
--- src/chkgmx.f	2016-05-09 12:15:45 +0000
+++ src/chkgmx.F90	2018-02-19 19:59:41 +0000
@@ -1,116 +1,104 @@
 !!@LICENSE
 
-      MODULE m_chkgmx
-
-! Used module procedures
-      use sys,       only: die     ! Termination routine
-      use m_minvec,  only: minvec  ! Finds a lattice basis of minimal length
-      use cellsubs,  only: reclat  ! Finds reciprocal unit cell
-
-! Used module parameters and variables
-      use precision, only: dp      ! Double precision real kind
-!      use parallel,  only: Node    ! My processor index
-
-      implicit none
-
-! Public procedures provided:
-      PUBLIC::
-     .  chkgmx,  ! Checks that a given cutoff is consistent with a mesh
-     .  meshKcut ! Returns the planewave cutoff of a periodic mesh
-
-! Public parameters, variables, and arrays:
-!     none
-
-      PRIVATE ! Nothing is declared public beyond this point
-
-      CONTAINS
-!---------------------------------------------------------------------------
-
-      SUBROUTINE CHKGMX (K,BG,MESH,G2MAX)
-
-      implicit none
-      real(dp),intent(in) :: K(3)    ! Vector in reciprocal unit cell (BZ)
-      real(dp),intent(in) :: BG(3,3) ! Reciprocal lattice vectors BG(ixyz,ivec)
-      integer, intent(in) :: MESH(3) ! Number of mesh point divisions of each 
-                                     ! real-space lattice vector.
-      real(dp),intent(inout):: G2MAX ! Planewave cutoff, which is reduced to
-                                     ! the maximum value supported by the mesh
-                                     ! (without aliasing), if this is smaller.
-
-      real(dp), parameter :: ZERO=0.0_dp,HALF=0.5_dp,
-     .                       TINY=1.0e-8_dp,BIG=1.0e20_dp
-
-      integer :: i, j, i1, i2, i3
-      real(dp):: bm(3,3), baux(3,3), ctransf(3,3),
-     .           g(3), gmod, gmax, g2msh, r
-
-      DO I=1,3
-        DO J=1,3
-          BM(J,I)=BG(J,I)*MESH(I)
-          Baux(J,I)=BG(J,I)*MESH(I)
-        ENDDO
-      ENDDO
-!
-!     Use Baux to avoid aliasing of arguments, as one is intent(in)
-!     and the other intent(out)...
-!
-      CALL MINVEC (Baux,BM, ctransf)
-      GMAX=BIG
-      DO I3=-1,1
-        DO I2=-1,1
-          DO I1=-1,1
-            IF (I1.EQ.0.AND.I2.EQ.0.AND.I3.EQ.0) GO TO 20
-              G(1)=BM(1,1)*I1+BM(1,2)*I2+BM(1,3)*I3
-              G(2)=BM(2,1)*I1+BM(2,2)*I2+BM(2,3)*I3
-              G(3)=BM(3,1)*I1+BM(3,2)*I2+BM(3,3)*I3
-              GMOD=SQRT(SUM(G**2))
-              R=HALF*GMOD-SUM(K*G)/GMOD
-              GMAX=MIN(GMAX,R)
-  20        CONTINUE
-          ENDDO
-        ENDDO
-      ENDDO
-      IF (GMAX.LT.ZERO) call die('CHKGMX: K NOT IN FIRST BZ')
-      G2MSH=GMAX*GMAX-TINY
-      IF (G2MSH.LT.G2MAX) THEN
-!       if (Node.eq.0) then
-!         WRITE(6,*) 'CHKGMX: MESH TOO SPARSE. GMAX HAS BEEN REDUCED'
-!         WRITE(6,*) 'CHKGMX: OLD G2MAX =',G2MAX
-!         WRITE(6,*) 'CHKGMX: NEW G2MAX =',G2MSH
-!       endif
-        G2MAX=G2MSH
-      ENDIF
-      END SUBROUTINE CHKGMX
-
-!----------------------------------------------------------------------
-
-      FUNCTION meshKcut( cell, nMesh )
-
-      ! Finds mesh planewave cutoff
-
-      use precision, only: dp
-
-      implicit none
-      real(dp),intent(in):: cell(3,3) ! Unit cell vectors cell(iXYZ,iVector)
-      integer, intent(in):: nMesh(3)  ! Mesh divisions of each vector
-      real(dp)           :: meshKcut  ! Mesh wavevector cutoff
-
-      real(dp):: g2max, k0(3), kcell(3,3)
-
-      ! Find reciprocal cell vectors
-      call reclat( cell, kcell, 1 )
-
-      ! Initialize arguments to call chkgmx
-      k0(1:3) = 0._dp
-      g2max = huge(g2max)
-
-      ! chkgmx will reduce g2max to square of mesh cutoff
-      call chkgmx( k0, kcell, nMesh, g2max )
-
-      ! Return wavevector cutoff
-      meshKcut = sqrt(g2max)
-
-      END FUNCTION meshKcut
-
-      END MODULE m_chkgmx
+module m_chkgmx
+
+  ! Used module procedures
+  use sys,       only: die     ! Termination routine
+  use m_minvec,  only: minvec  ! Finds a lattice basis of minimal length
+  use cellsubs,  only: reclat  ! Finds reciprocal unit cell
+
+  ! Used module parameters and variables
+  use precision, only: dp      ! Double precision real kind
+
+  implicit none
+
+  private ! Nothing is declared public beyond this point
+
+  ! Public procedures provided:
+  public :: chkgmx  ! Checks that a given cutoff is consistent with a mesh
+  public :: meshKcut ! Returns the planewave cutoff of a periodic mesh
+
+contains
+
+  subroutine chkgmx(K,BG,MESH,G2MAX)
+
+    real(dp),intent(in) :: K(3)    ! Vector in reciprocal unit cell (BZ)
+    real(dp),intent(in) :: BG(3,3) ! Reciprocal lattice vectors BG(ixyz,ivec)
+    integer, intent(in) :: MESH(3) ! Number of mesh point divisions of each 
+    ! real-space lattice vector.
+    real(dp),intent(inout):: G2MAX ! Planewave cutoff, which is reduced to
+    ! the maximum value supported by the mesh
+    ! (without aliasing), if this is smaller.
+
+    real(dp), parameter :: ZERO=0.0_dp,HALF=0.5_dp, &
+        TINY=1.0e-8_dp, BIG=1.0e20_dp
+
+    integer :: i, j, i1, i2, i3
+    real(dp):: bm(3,3), baux(3,3), ctransf(3,3), &
+        g(3), gmod, gmax, g2msh, r
+    
+    DO I=1,3
+      DO J=1,3
+        BM(J,I)=BG(J,I)*MESH(I)
+        Baux(J,I)=BG(J,I)*MESH(I)
+      END DO
+    END DO
+    !
+    !     Use Baux to avoid aliasing of arguments, as one is intent(in)
+    !     and the other intent(out)...
+    !
+    CALL MINVEC (Baux,BM, ctransf)
+    GMAX=BIG
+    DO I3=-1,1
+      DO I2=-1,1
+        DO I1=-1,1
+          IF (I1.EQ.0.AND.I2.EQ.0.AND.I3.EQ.0) cycle
+          G(1)=BM(1,1)*I1+BM(1,2)*I2+BM(1,3)*I3
+          G(2)=BM(2,1)*I1+BM(2,2)*I2+BM(2,3)*I3
+          G(3)=BM(3,1)*I1+BM(3,2)*I2+BM(3,3)*I3
+          GMOD=SQRT(SUM(G**2))
+          R=HALF*GMOD-SUM(K*G)/GMOD
+          GMAX=MIN(GMAX,R)
+        END DO
+      END DO
+    END DO
+    IF (GMAX.LT.ZERO) call die('CHKGMX: K NOT IN FIRST BZ')
+    G2MSH=GMAX*GMAX-TINY
+    IF (G2MSH.LT.G2MAX) THEN
+      !       if (Node.eq.0) then
+      !         WRITE(6,*) 'CHKGMX: MESH TOO SPARSE. GMAX HAS BEEN REDUCED'
+      !         WRITE(6,*) 'CHKGMX: OLD G2MAX =',G2MAX
+      !         WRITE(6,*) 'CHKGMX: NEW G2MAX =',G2MSH
+      !       endif
+      G2MAX=G2MSH
+    ENDIF
+  end subroutine chkgmx
+
+  !----------------------------------------------------------------------
+
+  function meshKcut( cell, nMesh )
+
+    ! Finds mesh planewave cutoff
+    real(dp),intent(in) :: cell(3,3) ! Unit cell vectors cell(iXYZ,iVector)
+    integer, intent(in) :: nMesh(3)  ! Mesh divisions of each vector
+    real(dp)           :: meshKcut  ! Mesh wavevector cutoff
+
+    real(dp) :: g2max, k0(3), kcell(3,3)
+
+    ! Find reciprocal cell vectors
+    call reclat( cell, kcell, 1 )
+
+    ! Initialize arguments to call chkgmx
+    k0(1:3) = 0._dp
+    g2max = huge(g2max)
+
+    ! chkgmx will reduce g2max to square of mesh cutoff
+    call chkgmx( k0, kcell, nMesh, g2max )
+
+    ! Return wavevector cutoff
+    meshKcut = sqrt(g2max)
+
+  end function meshKcut
+
+end module m_chkgmx
 

=== renamed file 'src/ggaxc.F' => 'src/ggaxc.F90'
--- src/ggaxc.F	2017-10-30 11:35:11 +0000
+++ src/ggaxc.F90	2018-02-19 19:59:41 +0000
@@ -44,156 +44,151 @@
 !
 !******************************************************************************
 
-      MODULE m_ggaxc
-
-      ! Used module procedures
-      use sys,     only: die     ! Termination routine
-      USE m_ldaxc, only: exchng  ! Local exchange
-      USE m_ldaxc, only: pw92c   ! Perdew & Wang, PRB, 45, 13244 (1992) correl
-
-
-      ! Used module parameters
-      use precision, only : dp   ! Double precision real kind
+MODULE m_ggaxc
+
+  ! Used module procedures
+  use sys,     only: die     ! Termination routine
+  USE m_ldaxc, only: exchng  ! Local exchange
+  USE m_ldaxc, only: pw92c   ! Perdew & Wang, PRB, 45, 13244 (1992) correl
+
+  ! Used module parameters
+  use precision, only : dp   ! Double precision real kind
 
 #ifdef LIBXC
-      use xc_f90_types_m
-      use xc_f90_lib_m
+  use xc_f90_types_m
+  use xc_f90_lib_m
 #endif
 
-      implicit none
-
-      PUBLIC::  
-     .  ggaxc,        ! General subroutine for all coded GGA XC functionals
-     .  blypxc,       ! Becke-Lee-Yang-Parr (see subroutine blypxc)
-     .  pbexc,        ! Perdew, Burke & Ernzerhof, PRL 77, 3865 (1996)
-     .  pbesolxc,     ! Perdew et al, PRL, 100, 136406 (2008)
-     .  pw91xc,       ! Perdew & Wang, JCP, 100, 1290 (1994)
-     .  revpbexc,     ! GGA Zhang & Yang, PRL 80,890(1998)
-     .  rpbexc,       ! Hammer, Hansen & Norskov, PRB 59, 7413 (1999)
-     .  am05xc,       ! Mattsson & Armiento, PRB, 79, 155101 (2009)
-     .  wcxc,         ! Wu-Cohen (see subroutine wcxc)
-     .  pbeJsJrLOxc,  ! Reparametrizations of the PBE functional by
-     .  pbeJsJrHEGxc, !   L.S.Pedroza et al, PRB 79, 201106 (2009) and
-     .  pbeGcGxLOxc,  !   M.M.Odashima et al, JCTC 5, 798 (2009)
-     .  pbeGcGxHEGxc, ! using 4 different combinations of criteria
-     .  pw86x,        ! Perdew & Wang, PRB 33, 8800 (1986) (exchange only)
-     .  pw86rx,       ! pw86x refit: Murray, Lee & Langreth, JCTC 5, 2754 (2009)
-     .  b88x,         ! Becke, PRA 38, 3098 (1988) (exchange only)
-     .  b88kbmx,      ! Klimes et al, JPCM 22, 022201 (2009) (exchange only)
-     .  c09x,         ! Cooper, PRB 81, 161104 (2010) (exchange only)
-     .  bhx           ! Berland & Hyldgaard, PRB 89, 035412 (2014) (exch. only)
-
-      PRIVATE  ! Nothing is declared public beyond this point
-
-      CONTAINS
-
-      SUBROUTINE GGAXC( AUTHOR, IREL, nSpin, D, GD,
-     .                  EPSX, EPSC, dEXdD, dECdD, dEXdGD, dECdGD
+  implicit none
+
+  private  ! everything is private, unless otherwise declared
+
+  public :: ggaxc        ! General subroutine for all coded GGA XC functionals
+  public :: blypxc       ! Becke-Lee-Yang-Parr (see subroutine blypxc)
+  public :: pbexc        ! Perdew, Burke & Ernzerhof, PRL 77, 3865 (1996)
+  public :: pbesolxc     ! Perdew et al, PRL, 100, 136406 (2008)
+  public :: pw91xc       ! Perdew & Wang, JCP, 100, 1290 (1994)
+  public :: revpbexc     ! GGA Zhang & Yang, PRL 80,890(1998)
+  public :: rpbexc       ! Hammer, Hansen & Norskov, PRB 59, 7413 (1999)
+  public :: am05xc       ! Mattsson & Armiento, PRB, 79, 155101 (2009)
+  public :: wcxc         ! Wu-Cohen (see subroutine wcxc)
+  public :: pbeJsJrLOxc  ! Reparametrizations of the PBE functional by
+  public :: pbeJsJrHEGxc !   L.S.Pedroza et al, PRB 79, 201106 (2009) and
+  public :: pbeGcGxLOxc  !   M.M.Odashima et al, JCTC 5, 798 (2009)
+  public :: pbeGcGxHEGxc ! using 4 different combinations of criteria
+  public :: pw86x        ! Perdew & Wang, PRB 33, 8800 (1986) (exchange only)
+  public :: pw86rx       ! pw86x refit: Murray, Lee & Langreth, JCTC 5, 2754 (2009)
+  public :: b88x         ! Becke, PRA 38, 3098 (1988) (exchange only)
+  public :: b88kbmx      ! Klimes et al, JPCM 22, 022201 (2009) (exchange only)
+  public :: c09x         ! Cooper, PRB 81, 161104 (2010) (exchange only)
+  public :: bhx          ! Berland & Hyldgaard, PRB 89, 035412 (2014) (exch. only)
+
+contains
+
+  subroutine GGAXC( AUTHOR, IREL, nSpin, D, GD, &
+      EPSX, EPSC, dEXdD, dECdD, dEXdGD, dECdGD &
 #ifdef LIBXC
-     .              , is_libxc_in, xc_func, xc_info )
+      , is_libxc_in, xc_func, xc_info )
 #else
-     .               )
-#endif
-
-C Finds the exchange and correlation energies at a point, and their
-C derivatives with respect to density and density gradient, in the
-C Generalized Gradient Correction approximation.
-C Lengths in Bohr, energies in Hartrees
-C Written by L.C.Balbas and J.M.Soler, Dec'96.
-C Modified by V.M.Garcia-Suarez to include non-collinear spin. June 2002
-C Non collinear part rewritten by J.M.Soler. Sept. 2009
-C Interface with LibXC added by Micael Oliveira, Jul.2014
-
-      implicit none
-
-      ! Input
-      character(len=*),intent(in):: AUTHOR  ! GGA flavour (author initials)
-      integer, intent(in) :: IREL           ! Relativistic exchange? 0=no, 1=yes
-      integer, intent(in) :: nSpin          ! Number of spin components
-      real(dp),intent(in) :: D(nSpin)       ! Local electron (spin) density
-      real(dp),intent(in) :: GD(3,nSpin)    ! Gradient of electron density
-
-      ! Output
-      real(dp),intent(out):: EPSX           ! Exchange energy per electron
-      real(dp),intent(out):: EPSC           ! Correlation energy per electron
-      real(dp),intent(out):: dEXdD(nSpin)   ! dEx/dDens, Ex=Int(dens*epsX)
-      real(dp),intent(out):: dECdD(nSpin)   ! dEc/dDens
-      real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens)
-      real(dp),intent(out):: dECdGD(3,nSpin) ! dEc/dGrad(Dens)
-#ifdef LIBXC
-      logical, intent(in), optional    :: is_libxc_in
-      type(xc_f90_pointer_t), optional :: xc_func, xc_info
-#endif
-
-      ! Internal variables and arrays
-      integer    :: NS, is, ix
-      real(dp)   :: DD(2), dECdDD(2), dEXdDD(2),
-     .              dDDdD(2,4), dDTOTdD(4), dDPOLdD(4),
-     .              dECdGDD(3,2), dEXdGDD(3,2),
-     .              dGDDdD(3,2,4), dGDDdGD(2,4),
-     .              dGDTOTdD(3,4), dGDPOLdD(3,4),
-     .              dGDTOTdGD(4), dGDPOLdGD(4),
-     .              DPOL, DTOT, GDD(3,2), GDTOT(3), GDPOL(3),
-     .              VPOL
-      real(dp), parameter :: DENMIN = 1.e-12_dp
-
-      real(dp):: theta, phi, c2, s2, st, ct, cp, sp, dpolz, dpolxy
-
-      logical , parameter :: old_scheme = .true.
-
-#ifdef LIBXC
-      logical  :: is_libxc
-      real(dp) :: eps, dedn(nspin), sgm(3), dedsgm(3), dedgd(3, nspin)
-#endif
-
-      ! Handle non-collinear spin case
-      if (nSpin==4) then
-        NS = 2             ! Diagonal spin components
-
-        if ( old_scheme ) then
-         DTOT = D(1) + D(2)
-         dpolz= D(1) - D(2)
-         dpolxy= 2.0d0*sqrt(d(3)**2+d(4)**2)
-         dpol  = sqrt( dpolz**2 + dpolxy**2 )
-         if ( dpol.gt.1.0d-12 ) then
+    )
+#endif
+
+    ! Finds the exchange and correlation energies at a point, and their
+    ! derivatives with respect to density and density gradient, in the
+    ! Generalized Gradient Correction approximation.
+    ! Lengths in Bohr, energies in Hartrees
+    ! Written by L.C.Balbas and J.M.Soler, Dec'96.
+    ! Modified by V.M.Garcia-Suarez to include non-collinear spin. June 2002
+    ! Non collinear part rewritten by J.M.Soler. Sept. 2009
+    ! Interface with LibXC added by Micael Oliveira, Jul.2014
+
+    ! Input
+    character(len=*),intent(in):: AUTHOR  ! GGA flavour (author initials)
+    integer, intent(in) :: IREL           ! Relativistic exchange? 0=no, 1=yes
+    integer, intent(in) :: nSpin          ! Number of spin components
+    real(dp),intent(in) :: D(nSpin)       ! Local electron (spin) density
+    real(dp),intent(in) :: GD(3,nSpin)    ! Gradient of electron density
+
+    ! Output
+    real(dp),intent(out):: EPSX           ! Exchange energy per electron
+    real(dp),intent(out):: EPSC           ! Correlation energy per electron
+    real(dp),intent(out):: dEXdD(nSpin)   ! dEx/dDens, Ex=Int(dens*epsX)
+    real(dp),intent(out):: dECdD(nSpin)   ! dEc/dDens
+    real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens)
+    real(dp),intent(out):: dECdGD(3,nSpin) ! dEc/dGrad(Dens)
+#ifdef LIBXC
+    logical, intent(in), optional    :: is_libxc_in
+    type(xc_f90_pointer_t), optional :: xc_func, xc_info
+#endif
+
+    ! Internal variables and arrays
+    integer    :: NS, is, ix
+    real(dp)   :: DD(2), dECdDD(2), dEXdDD(2), &
+        dDDdD(2,4), dDTOTdD(4), dDPOLdD(4), &
+        dECdGDD(3,2), dEXdGDD(3,2), &
+        dGDDdD(3,2,4), dGDDdGD(2,4), &
+        dGDTOTdD(3,4), dGDPOLdD(3,4), &
+        dGDTOTdGD(4), dGDPOLdGD(4), &
+        DPOL, DTOT, GDD(3,2), GDTOT(3), GDPOL(3), &
+        VPOL
+    real(dp), parameter :: DENMIN = 1.e-12_dp
+
+    real(dp):: theta, phi, c2, s2, st, ct, cp, sp, dpolz, dpolxy
+
+    logical , parameter :: old_scheme = .true.
+
+#ifdef LIBXC
+    logical  :: is_libxc
+    real(dp) :: eps, dedn(nspin), sgm(3), dedsgm(3), dedgd(3, nspin)
+#endif
+
+    ! Handle non-collinear spin case
+    if (nSpin==4) then
+      NS = 2             ! Diagonal spin components
+
+      if ( old_scheme ) then
+        DTOT = D(1) + D(2)
+        dpolz= D(1) - D(2)
+        dpolxy= 2.0d0*sqrt(d(3)**2+d(4)**2)
+        dpol  = sqrt( dpolz**2 + dpolxy**2 )
+        if ( dpol.gt.1.0d-12 ) then
           theta = atan2(dpolxy,dpolz)
-         else
+        else
           theta = 0.0_dp
-         endif
-         C2 = COS(THETA/2.0_dp)
-         S2 = SIN(THETA/2.0_dp)
-         ST = SIN(THETA)
-         CT = COS(THETA)
-         PHI = ATAN2(-D(4),D(3))
-         CP = COS(PHI)
-         SP = SIN(PHI)
-
-         DD(1) = 0.5D0 * ( DTOT + DPOL )
-         DD(2) = 0.5D0 * ( DTOT - DPOL )
-         DO IX = 1,3
-          GDD(IX,1) = GD(IX,1)*C2**2 + GD(IX,2)*S2**2 +
-     .                2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP)
-          GDD(IX,2) = GD(IX,1)*S2**2 + GD(IX,2)*C2**2 -
-     .                2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP)
-         ENDDO
-
-        else
+        endif
+        C2 = COS(THETA/2.0_dp)
+        S2 = SIN(THETA/2.0_dp)
+        ST = SIN(THETA)
+        CT = COS(THETA)
+        PHI = ATAN2(-D(4),D(3))
+        CP = COS(PHI)
+        SP = SIN(PHI)
+
+        DD(1) = 0.5D0 * ( DTOT + DPOL )
+        DD(2) = 0.5D0 * ( DTOT - DPOL )
+        DO IX = 1,3
+          GDD(IX,1) = GD(IX,1)*C2**2 + GD(IX,2)*S2**2 + &
+              2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP)
+          GDD(IX,2) = GD(IX,1)*S2**2 + GD(IX,2)*C2**2 - &
+              2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP)
+        ENDDO
+
+      else
 
         ! Find eigenvalues of density matrix Dij (diagonal densities DD, i.e.
         ! up and down densities along the spin direction). Note convention: 
         ! D(1)=D11, D(2)=D22, D(3)=Re(D12)=Re(D21), D(4)=Im(D12)=-Im(D21)
         DTOT = D(1) + D(2)                           ! DensTot (DensUp+DensDn)
-        DPOL = SQRT( (D(1)-D(2))**2                  ! DensPol (DensUp-DensDn)
-     .              + 4*(D(3)**2+D(4)**2) )
+        DPOL = SQRT( (D(1)-D(2))**2 &                ! DensPol (DensUp-DensDn)
+            + 4*(D(3)**2+D(4)**2) )
         DD(1) = ( DTOT + DPOL ) / 2                  ! DensUp
         DD(2) = ( DTOT - DPOL ) / 2                  ! DensDn
         DPOL = max( DPOL , DENMIN )                  ! Avoid division by zero
 
         ! Find gradients of up and down densities
         GDTOT(:) = GD(:,1) + GD(:,2)                 ! Grad(DensTot)
-        GDPOL(:) = ( (D(1)-D(2))*(GD(:,1)-GD(:,2))   ! Grad(DensPol)
-     .             + 4*(D(3)*GD(:,3)+D(4)*GD(:,4)) )
-     .             / DPOL
+        GDPOL(:) = ( (D(1)-D(2))*(GD(:,1)-GD(:,2)) & ! Grad(DensPol)
+            + 4*(D(3)*GD(:,3)+D(4)*GD(:,4)) ) / DPOL
         GDD(:,1) = ( GDTOT(:) + GDPOL(:) ) / 2       ! Grad(DensUp)
         GDD(:,2) = ( GDTOT(:) - GDPOL(:) ) / 2       ! Grad(DensDn)
 
@@ -209,18 +204,18 @@
 
         ! Derivatives of grad(Dup) and grad(Ddn) with respect to D(i)
         dGDTOTdD(1:3,1:4) = 0                        ! dGradDensTot/dD(i)
-        dGDPOLdD(:,1) = + (GD(:,1)-GD(:,2))/DPOL     ! dGradDensPol/dD(i)
-     .                  - GDPOL(:) * dDPOLdD(1)/DPOL
-        dGDPOLdD(:,2) = - (GD(:,1)-GD(:,2))/DPOL
-     .                  - GDPOL(:) * dDPOLdD(2)/DPOL
-        dGDPOLdD(:,3) = 4*GD(:,3)/DPOL
-     .                  - GDPOL(:) * dDPOLdD(3)/DPOL
-        dGDPOLdD(:,4) = 4*GD(:,4)/DPOL
-     .                  - GDPOL(:) * dDPOLdD(4)/DPOL
-        dGDDdD(:,1,:) = ( dGDTOTdD(:,:)              ! dGradDensUp/dD(i)
-     .                  + dGDPOLdD(:,:) ) / 2
-        dGDDdD(:,2,:) = ( dGDTOTdD(:,:)              ! dGradDensDn/dD(i)
-     .                  - dGDPOLdD(:,:) ) / 2
+        dGDPOLdD(:,1) = + (GD(:,1)-GD(:,2))/DPOL &   ! dGradDensPol/dD(i)
+            - GDPOL(:) * dDPOLdD(1)/DPOL
+        dGDPOLdD(:,2) = - (GD(:,1)-GD(:,2))/DPOL &
+            - GDPOL(:) * dDPOLdD(2)/DPOL
+        dGDPOLdD(:,3) = 4*GD(:,3)/DPOL &
+            - GDPOL(:) * dDPOLdD(3)/DPOL
+        dGDPOLdD(:,4) = 4*GD(:,4)/DPOL &
+            - GDPOL(:) * dDPOLdD(4)/DPOL
+        dGDDdD(:,1,:) = ( dGDTOTdD(:,:) &             ! dGradDensUp/dD(i)
+            + dGDPOLdD(:,:) ) / 2
+        dGDDdD(:,2,:) = ( dGDTOTdD(:,:) &             ! dGradDensDn/dD(i)
+            - dGDPOLdD(:,:) ) / 2
 
         ! Derivatives of grad(Dup) and grad(Ddn) with respect to grad(D(i))
         dGDTOTdGD(1:2) = 1                           ! dGradDensTot/dGradD(i)
@@ -229,190 +224,188 @@
         dGDPOLdGD(2) = -( D(1) - D(2) ) / DPOL
         dGDPOLdGD(3) = 4 * D(3) / DPOL
         dGDPOLdGD(4) = 4 * D(4) / DPOL
-        dGDDdGD(1,:) = ( dGDTOTdGD(:)                ! dGradDensUp/dGradD(i)
-     .                 + dGDPOLdGD(:) ) / 2
-        dGDDdGD(2,:) = ( dGDTOTdGD(:)                ! dGradDensDn/dGradD(i)
-     .                 - dGDPOLdGD(:) ) / 2
-       endif
-        
-      else if (nSpin==1 .or. nSpin==2) then ! Normal (collinear) spin
-        NS = nSpin
-        DD(1:NS) = max( D(1:NS), 0.0_dp ) ! ag: Avoid negative densities
-        GDD(1:3,1:NS) = GD(1:3,1:NS)
-      else
-        call die('ggaxc: ERROR: invalid value of nSpin')
-      end if ! (nSpin==4)
-
-      ! Select functional to find energy density and its derivatives
+        dGDDdGD(1,:) = ( dGDTOTdGD(:) &              ! dGradDensUp/dGradD(i)
+            + dGDPOLdGD(:) ) / 2
+        dGDDdGD(2,:) = ( dGDTOTdGD(:) &              ! dGradDensDn/dGradD(i)
+            - dGDPOLdGD(:) ) / 2
+      endif
+
+    else if (nSpin==1 .or. nSpin==2) then ! Normal (collinear) spin
+      NS = nSpin
+      DD(1:NS) = max( D(1:NS), 0.0_dp ) ! ag: Avoid negative densities
+      GDD(1:3,1:NS) = GD(1:3,1:NS)
+    else
+      call die('ggaxc: ERROR: invalid value of nSpin')
+    end if ! (nSpin==4)
+
+    ! Select functional to find energy density and its derivatives
 #ifdef LIBXC
-      is_libxc = .false.
-      if (present(is_libxc_in)) then
-         if (is_libxc_in) then
-            if ((.not. present(xc_func)) .or.
-     $           (.not. present(xc_info))) then
-               call die("xc_func and xc_info not present")
-            endif
-            is_libxc = .true.
-         endif
+    is_libxc = .false.
+    if (present(is_libxc_in)) then
+      if (is_libxc_in) then
+        if ((.not. present(xc_func)) .or. &
+            (.not. present(xc_info))) then
+          call die("xc_func and xc_info not present")
+        endif
+        is_libxc = .true.
       endif
-
-      if (is_libxc) then
-        sgm(1) = sum(GDD(1:3,1)*GDD(1:3,1))
-        IF (nspin == 1) THEN
-           ! do nothing
-        ELSE
-          sgm(2) = sum(GDD(1:3,1)*GDD(1:3,2)) 
-          sgm(3) = sum(GDD(1:3,2)*GDD(1:3,2)) 
-        ENDIF
-        IF (xc_f90_info_family(xc_info) /= XC_FAMILY_GGA) THEN
-          call die('GGAXC: Functional is not a GGA')
-        ENDIF
-
-        call xc_f90_gga_exc_vxc(xc_func, 1, DD(1), sgm(1), eps, dedn(1), 
-     .                          dedsgm(1))
-        IF (nspin == 1) THEN
-          dedgd(1:3, 1) = 2.0D0*dedsgm(1)*GDD(1:3, 1)
-        ELSE
-          dedgd(1:3, 1) = 2.0D0*dedsgm(1)*GDD(1:3, 1) + 
-     .                    dedsgm(2)*GDD(1:3, 2)
-          dedgd(1:3, 2) = 2.0D0*dedsgm(3)*GDD(1:3, 2) + 
-     .                    dedsgm(2)*GDD(1:3, 1)
-        ENDIF
-        IF (xc_f90_info_kind(xc_info) == XC_CORRELATION) THEN
-          EPSC = eps
-          dECdDD(1:nspin) = dedn(1:nspin)
-          dECdGDD(1:3, 1:nspin) = dedgd(1:3, 1:nspin)
-          EPSX = 0.0_dp
-          dEXdDD(1:nspin) = 0.0_dp
-          dEXdGDD(1:3, 1:nspin) = 0.0_dp
-        ELSE if (xc_f90_info_kind(xc_info) == XC_EXCHANGE) THEN
-          EPSX = eps
-          dEXdDD(1:nspin) = dedn(1:nspin)
-          dEXdGDD(1:3, 1:nspin) = dedgd(1:3, 1:nspin)
-          EPSC = 0.0_dp
-          dECdDD(1:nspin) = 0.0_dp
-          dECdGDD(1:3, 1:nspin) = 0.0_dp
-        ELSE ! combined functional, use an arbitrary 50/50 split
-          EPSX = 0.5_dp * eps
-          dEXdDD(1:nspin) = 0.5_dp * dedn(1:nspin)
-          dEXdGDD(1:3, 1:nspin) = 0.5_dp * dedgd(1:3, 1:nspin)
-          !
-          EPSC = EPSX
-          dECdDD(1:nspin) = dEXdDD(1:nspin)
-          dECdGDD(1:3, 1:nspin) = dEXdGDD(1:3,1:nspin)
-         ENDIF
-
-      else IF (AUTHOR.EQ.'PBE' .OR. AUTHOR.EQ.'pbe') THEN
+    endif
+
+    if (is_libxc) then
+      sgm(1) = sum(GDD(1:3,1)*GDD(1:3,1))
+      IF (nspin == 1) THEN
+        ! do nothing
+      ELSE
+        sgm(2) = sum(GDD(1:3,1)*GDD(1:3,2)) 
+        sgm(3) = sum(GDD(1:3,2)*GDD(1:3,2)) 
+      ENDIF
+      IF (xc_f90_info_family(xc_info) /= XC_FAMILY_GGA) THEN
+        call die('GGAXC: Functional is not a GGA')
+      ENDIF
+
+      call xc_f90_gga_exc_vxc(xc_func, 1, DD(1), sgm(1), eps, dedn(1), &
+          dedsgm(1))
+      IF (nspin == 1) THEN
+        dedgd(1:3, 1) = 2.0D0*dedsgm(1)*GDD(1:3, 1)
+      ELSE
+        dedgd(1:3, 1) = 2.0D0*dedsgm(1)*GDD(1:3, 1) + dedsgm(2)*GDD(1:3, 2)
+        dedgd(1:3, 2) = 2.0D0*dedsgm(3)*GDD(1:3, 2) + dedsgm(2)*GDD(1:3, 1)
+      ENDIF
+      IF (xc_f90_info_kind(xc_info) == XC_CORRELATION) THEN
+        EPSC = eps
+        dECdDD(1:nspin) = dedn(1:nspin)
+        dECdGDD(1:3, 1:nspin) = dedgd(1:3, 1:nspin)
+        EPSX = 0.0_dp
+        dEXdDD(1:nspin) = 0.0_dp
+        dEXdGDD(1:3, 1:nspin) = 0.0_dp
+      ELSE if (xc_f90_info_kind(xc_info) == XC_EXCHANGE) THEN
+        EPSX = eps
+        dEXdDD(1:nspin) = dedn(1:nspin)
+        dEXdGDD(1:3, 1:nspin) = dedgd(1:3, 1:nspin)
+        EPSC = 0.0_dp
+        dECdDD(1:nspin) = 0.0_dp
+        dECdGDD(1:3, 1:nspin) = 0.0_dp
+      ELSE ! combined functional, use an arbitrary 50/50 split
+        EPSX = 0.5_dp * eps
+        dEXdDD(1:nspin) = 0.5_dp * dedn(1:nspin)
+        dEXdGDD(1:3, 1:nspin) = 0.5_dp * dedgd(1:3, 1:nspin)
+        !
+        EPSC = EPSX
+        dECdDD(1:nspin) = dEXdDD(1:nspin)
+        dECdGDD(1:3, 1:nspin) = dEXdGDD(1:3,1:nspin)
+      ENDIF
+
+    else IF (AUTHOR.EQ.'PBE' .OR. AUTHOR.EQ.'pbe') THEN
 #else         
-      IF (AUTHOR.EQ.'PBE' .OR. AUTHOR.EQ.'pbe') THEN
+    if (AUTHOR.EQ.'PBE' .OR. AUTHOR.EQ.'pbe') THEN
 #endif         
-        CALL PBEXC( IREL, NS, DD, GDD,                  ! JMS
-     .              EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSE IF (AUTHOR.EQ.'RPBE'.OR.AUTHOR.EQ.'rpbe') THEN
-        CALL RPBEXC( IREL, NS, DD, GDD,                 ! MVFS
-     .               EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSE IF (AUTHOR.EQ.'WC'.OR.AUTHOR.EQ.'wc') THEN
-        CALL WCXC( IREL, NS, DD, GDD,                   ! MVFS
-     .               EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSE IF (AUTHOR.EQ.'REVPBE'.OR.AUTHOR.EQ.'revpbe'
-     .                           .OR.AUTHOR.EQ.'revPBE') THEN
-        CALL REVPBEXC( IREL, NS, DD, GDD,               ! EA
-     .               EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSE IF (AUTHOR.EQ.'BLYP'.OR.AUTHOR.EQ.'blyp'.OR.
-     .         AUTHOR.EQ.'LYP'.OR.AUTHOR.EQ.'lyp') THEN
-        CALL BLYPXC( NS, DD, GDD,                       ! AG
-     .               EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'PW91' .OR. AUTHOR.EQ.'pw91') THEN
-        CALL PW91XC( IREL, NS, DD, GDD,                 ! AG
-     .               EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'PBESOL' .OR. AUTHOR.EQ.'pbesol' .OR.
-     .        AUTHOR.EQ.'PBEsol') THEN
-        CALL PBESOLXC( IREL, NS, DD, GDD,
-     .                 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'AM05' .OR. AUTHOR.EQ.'am05') THEN
-        CALL AM05XC( IREL, NS, DD, GDD,
-     .               EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'PBEJsJrLO' .OR.
-     .        AUTHOR.EQ.'pbejsjrlo' .OR.
-     .        AUTHOR.EQ.'PBEJSJRLO') THEN
-        CALL PBEJsJrLOxc( IREL, NS, DD, GDD,
-     .                  EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'PBEJsJrHEG' .OR.
-     .        AUTHOR.EQ.'pbejsjrheg' .OR.
-     .        AUTHOR.EQ.'PBEJSJRHEG') THEN
-        CALL PBEJsJrHEGxc( IREL, NS, DD, GDD,
-     .                  EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'PBEGcGxLO' .OR.
-     .        AUTHOR.EQ.'pbegcgxlo' .OR.
-     .        AUTHOR.EQ.'PBEGCGXLO') THEN
-        CALL PBEGcGxLOxc( IREL, NS, DD, GDD,
-     .                  EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'PBEGcGxHEG' .OR.
-     .        AUTHOR.EQ.'pbegcgxheg' .OR.
-     .        AUTHOR.EQ.'PBEGCGXHEG') THEN
-        CALL PBEGcGxHEGxc( IREL, NS, DD, GDD,
-     .                  EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'PW86' .OR. AUTHOR.EQ.'pw86') THEN
-        CALL PW86X( IREL, NS, DD, GDD,
-     .              EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'PW86R' .OR. AUTHOR.EQ.'pw86r') THEN
-        CALL PW86RX( IREL, NS, DD, GDD,
-     .               EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'B88' .OR. AUTHOR.EQ.'b88') THEN
-        CALL B88X( IREL, NS, DD, GDD,
-     .             EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'B88KBM' .OR. AUTHOR.EQ.'b88kbm') THEN
-        CALL B88KBMX( IREL, NS, DD, GDD,
-     .             EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'C09' .OR. AUTHOR.EQ.'c09') THEN
-        CALL C09X( IREL, NS, DD, GDD,
-     .             EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
-
-      ELSEIF (AUTHOR.EQ.'BH' .OR. AUTHOR.EQ.'bh') THEN
-        CALL BHX( IREL, NS, DD, GDD,
-     .            EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 
-
-      ELSE
-        call die('GGAXC: Unknown author ' // trim(AUTHOR))
-      ENDIF
-
-      ! Find dE/dD(i) and dE/dGradD(i). Note convention:
-      ! DEDD(1)=dE/dD11, DEDD(2)=dE/dD22, DEDD(3)=Re(dE/dD12)=Re(dE/dD21),
-      ! DEDD(4)=Im(dE/dD12)=-Im(dE/D21)
-      if (nSpin==4) then  ! Non colinear spin
-        ! dE/dD(i) = dE/dDup * dDup/dD(i) + dE/dDdn * dDdn/dD(i)
-        !          + dE/dGDup * dGDup/dD(i) + dE/dGDdn * dGDdn/dD(i)
-        ! dE/dGradD(i) = dE/dGDup * dGDup/dGD(i) + dE/dGDdn * dGDdn/dGD(i)
-
-        if ( old_scheme ) then
-         VPOL  = (dExdDD(1)-dExdDD(2)) * ct
-         dExdD(1) = 0.5D0 * ( dExdDD(1) + dExdDD(2) + VPOL )
-         dExdD(2) = 0.5D0 * ( dExdDD(1) + dExdDD(2) - VPOL )
-         dExdD(3) = 0.5d0 * (dExdDD(1)-dExdDD(2)) * st * cp
-         dExdD(4) =-0.5d0 * (dExdDD(1)-dExdDD(2)) * st * sp
-         VPOL  = (dEcdDD(1)-dEcdDD(2)) * ct
-         dEcdD(1) = 0.5D0 * ( dEcdDD(1) + dEcdDD(2) + VPOL )
-         dEcdD(2) = 0.5D0 * ( dEcdDD(1) + dEcdDD(2) - VPOL )
-         dEcdD(3) = 0.5d0 * (dEcdDD(1)-dEcdDD(2)) * st * cp
-         dEcdD(4) =-0.5d0 * (dEcdDD(1)-dEcdDD(2)) * st * sp
-C Gradient terms
-         DO IX = 1,3
+      CALL PBEXC( IREL, NS, DD, GDD, &                 ! JMS
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSE IF (AUTHOR.EQ.'RPBE'.OR.AUTHOR.EQ.'rpbe') THEN
+      CALL RPBEXC( IREL, NS, DD, GDD, &                ! MVFS
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSE IF (AUTHOR.EQ.'WC'.OR.AUTHOR.EQ.'wc') THEN
+      CALL WCXC( IREL, NS, DD, GDD, &                  ! MVFS
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSE IF (AUTHOR.EQ.'REVPBE'.OR.AUTHOR.EQ.'revpbe' &
+        .OR.AUTHOR.EQ.'revPBE') THEN
+      CALL REVPBEXC( IREL, NS, DD, GDD, &              ! EA
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSE IF (AUTHOR.EQ.'BLYP'.OR.AUTHOR.EQ.'blyp'.OR. &
+        AUTHOR.EQ.'LYP'.OR.AUTHOR.EQ.'lyp') THEN
+      CALL BLYPXC( NS, DD, GDD, &                       ! AG
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'PW91' .OR. AUTHOR.EQ.'pw91') THEN
+      CALL PW91XC( IREL, NS, DD, GDD, &                ! AG
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'PBESOL' .OR. AUTHOR.EQ.'pbesol' .OR. &
+        AUTHOR.EQ.'PBEsol') THEN
+      CALL PBESOLXC( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'AM05' .OR. AUTHOR.EQ.'am05') THEN
+      CALL AM05XC( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'PBEJsJrLO' .OR. &
+        AUTHOR.EQ.'pbejsjrlo' .OR. &
+        AUTHOR.EQ.'PBEJSJRLO') THEN
+      CALL PBEJsJrLOxc( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'PBEJsJrHEG' .OR. &
+        AUTHOR.EQ.'pbejsjrheg' .OR. &
+        AUTHOR.EQ.'PBEJSJRHEG') THEN
+      CALL PBEJsJrHEGxc( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'PBEGcGxLO' .OR. &
+        AUTHOR.EQ.'pbegcgxlo' .OR. &
+        AUTHOR.EQ.'PBEGCGXLO') THEN
+      CALL PBEGcGxLOxc( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'PBEGcGxHEG' .OR. &
+        AUTHOR.EQ.'pbegcgxheg' .OR. &
+        AUTHOR.EQ.'PBEGCGXHEG') THEN
+      CALL PBEGcGxHEGxc( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'PW86' .OR. AUTHOR.EQ.'pw86') THEN
+      CALL PW86X( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'PW86R' .OR. AUTHOR.EQ.'pw86r') THEN
+      CALL PW86RX( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'B88' .OR. AUTHOR.EQ.'b88') THEN
+      CALL B88X( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'B88KBM' .OR. AUTHOR.EQ.'b88kbm') THEN
+      CALL B88KBMX( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'C09' .OR. AUTHOR.EQ.'c09') THEN
+      CALL C09X( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
+
+    ELSEIF (AUTHOR.EQ.'BH' .OR. AUTHOR.EQ.'bh') THEN
+      CALL BHX( IREL, NS, DD, GDD, &
+          EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 
+
+    ELSE
+      call die('GGAXC: Unknown author ' // trim(AUTHOR))
+    ENDIF
+
+    ! Find dE/dD(i) and dE/dGradD(i). Note convention:
+    ! DEDD(1)=dE/dD11, DEDD(2)=dE/dD22, DEDD(3)=Re(dE/dD12)=Re(dE/dD21),
+    ! DEDD(4)=Im(dE/dD12)=-Im(dE/D21)
+    if (nSpin==4) then  ! Non colinear spin
+      ! dE/dD(i) = dE/dDup * dDup/dD(i) + dE/dDdn * dDdn/dD(i)
+      !          + dE/dGDup * dGDup/dD(i) + dE/dGDdn * dGDdn/dD(i)
+      ! dE/dGradD(i) = dE/dGDup * dGDup/dGD(i) + dE/dGDdn * dGDdn/dGD(i)
+
+      if ( old_scheme ) then
+        VPOL  = (dExdDD(1)-dExdDD(2)) * ct
+        dExdD(1) = 0.5D0 * ( dExdDD(1) + dExdDD(2) + VPOL )
+        dExdD(2) = 0.5D0 * ( dExdDD(1) + dExdDD(2) - VPOL )
+        dExdD(3) = 0.5d0 * (dExdDD(1)-dExdDD(2)) * st * cp
+        dExdD(4) =-0.5d0 * (dExdDD(1)-dExdDD(2)) * st * sp
+        VPOL  = (dEcdDD(1)-dEcdDD(2)) * ct
+        dEcdD(1) = 0.5D0 * ( dEcdDD(1) + dEcdDD(2) + VPOL )
+        dEcdD(2) = 0.5D0 * ( dEcdDD(1) + dEcdDD(2) - VPOL )
+        dEcdD(3) = 0.5d0 * (dEcdDD(1)-dEcdDD(2)) * st * cp
+        dEcdD(4) =-0.5d0 * (dEcdDD(1)-dEcdDD(2)) * st * sp
+        ! Gradient terms
+        DO IX = 1,3
           dExdGD(IX,1) = dExdGDD(IX,1)*C2**2 + dExdGDD(IX,2)*S2**2
           dExdGD(IX,2) = dExdGDD(IX,1)*S2**2 + dExdGDD(IX,2)*C2**2
           dExdGD(IX,3) = 0.5D0*(dExdGDD(IX,1) - dExdGDD(IX,2))*ST*CP
@@ -421,15 +414,15 @@
           dEcdGD(IX,2) = dEcdGDD(IX,1)*S2**2 + dEcdGDD(IX,2)*C2**2
           dEcdGD(IX,3) = 0.5D0*(dEcdGDD(IX,1) - dEcdGDD(IX,2))*ST*CP
           dEcdGD(IX,4) =-0.5D0*(dEcdGDD(IX,1) - dEcdGDD(IX,2))*ST*SP
-         enddo
+        enddo
 
-        else
+      else
 
         do is = 1,4
-          dEXdD(is) = sum( dEXdDD(:) * dDDdD(:,is) )
-     .              + sum( dEXdGDD(:,:) * dGDDdD(:,:,is) )
-          dECdD(is) = sum( dECdDD(:) * dDDdD(:,is) )
-     .              + sum( dECdGDD(:,:) * dGDDdD(:,:,is) )
+          dEXdD(is) = sum( dEXdDD(:) * dDDdD(:,is) ) &
+              + sum( dEXdGDD(:,:) * dGDDdD(:,:,is) )
+          dECdD(is) = sum( dECdDD(:) * dDDdD(:,is) ) &
+              + sum( dECdGDD(:,:) * dGDDdD(:,:,is) )
           do ix = 1,3
             dEXdGD(ix,is) = sum( dEXdGDD(ix,:) * dGDDdGD(:,is) )
             dECdGD(ix,is) = sum( dECdGDD(ix,:) * dGDDdGD(:,is) )
@@ -445,881 +438,858 @@
         dEXdGD(:,3:4) = dEXdGD(:,3:4) / 2
         dECdGD(:,3:4) = dECdGD(:,3:4) / 2
 
-       endif
-      else   ! Collinear spin => just copy derivatives to output arrays
-        dEXdD(1:nSpin) = dEXdDD(1:nSpin)
-        dECdD(1:nSpin) = dECdDD(1:nSpin)
-        dEXdGD(:,1:nSpin) = dEXdGDD(:,1:nSpin)
-        dECdGD(:,1:nSpin) = dECdGDD(:,1:nSpin)
-      end if ! (nSpin==4)
-
-      END SUBROUTINE GGAXC
-
-
-
-      SUBROUTINE PBEformXC( beta, mu, kappa, iRel, nSpin, Dens, GDens,
-     .                      EX, EC, dEXdD, dECdD, dEXdGD, dECdGD )
-
-C *********************************************************************
-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
-C functional form, but with variable values for parameters beta, mu, and
-C kappa. Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
-C Written by L.C.Balbas and J.M.Soler. December 1996. 
-C ******** INPUT ******************************************************
-C REAL*8  BETA           : Parameter beta of the PBE functional
-C REAL*8  MU             : Parameter mu of the PBE functional
-C REAL*8  KAPPA          : Parameter kappa of the PBE functional
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C EXCHNG, PW92C
-C ********************************************************************
-
-      implicit none
-
-C Input
-      real(dp),intent(in) :: beta           ! Parameter of PBE functional
-      real(dp),intent(in) :: mu             ! Parameter of PBE functional
-      real(dp),intent(in) :: kappa          ! Parameter of PBE functional
-      integer, intent(in) :: iRel           ! Relativistic exchange? 0=no, 1=yes
-      integer, intent(in) :: nSpin          ! Number of spin components
-      real(dp),intent(in) :: Dens(nSpin)    ! Local electron (spin) density
-      real(dp),intent(in) :: GDens(3,nSpin) ! Gradient of electron density
-
-C Output
-      real(dp),intent(out):: EX             ! Exchange energy per electron
-      real(dp),intent(out):: EC             ! Correlation energy per electron
-      real(dp),intent(out):: dEXdD(nSpin)   ! dEx/dDens, Ex=Int(dens*epsX)
-      real(dp),intent(out):: dECdD(nSpin)   ! dEc/dDens
-      real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens)
-      real(dp),intent(out):: dECdGD(3,nSpin) ! dEc/dGrad(Dens)
-
-C Internal variables
-      INTEGER
-     .  IS, IX
-      real(dp)
-     .  A, D(2), DADD, DECUDD, DENMIN, 
-     .  DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
-     .  DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
-     .  DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, 
-     .  DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), 
-     .  ECUNIF, EXUNIF,
-     .  F, F1, F2, F3, F4, FC, FX, FOUTHD,
-     .  GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
-     .  H, HALF, KF, KFS, KS, PHI, PI, RS, S,
-     .  T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      PARAMETER ( DENMIN = 1.D-12 )
-      PARAMETER ( GDMIN  = 1.D-12 )
-
-C Fix some numerical parameters
-      PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
-     .            THD=1.D0/3.D0, THRHLF=1.5D0,
-     .            TWO=2.D0, TWOTHD=2.D0/3.D0 )
-
-C Fix some more numerical constants
-      PI = 4 * ATAN(1.D0)
-      GAMMA = (1 - LOG(TWO)) / PI**2
-
-C Translate density and its gradient to new variables
-      IF (nspin .EQ. 1) THEN
-        D(1) = HALF * Dens(1)
-        D(2) = D(1)
-        DT = MAX( DENMIN, Dens(1) )
-        DO 10 IX = 1,3
-          GD(IX,1) = HALF * GDens(IX,1)
-          GD(IX,2) = GD(IX,1)
-          GDT(IX) = GDens(IX,1)
-   10   CONTINUE
-      ELSE
-        D(1) = Dens(1)
-        D(2) = Dens(2)
-        DT = MAX( DENMIN, Dens(1)+Dens(2) )
-        DO 20 IX = 1,3
-          GD(IX,1) = GDens(IX,1)
-          GD(IX,2) = GDens(IX,2)
-          GDT(IX) = GDens(IX,1) + GDens(IX,2)
-   20   CONTINUE
-      ENDIF
-      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = MAX( GDMIN, GDMT )
-
-C Find local correlation energy and potential
-      CALL PW92C( 2, D, ECUNIF, VCUNIF )
-
-C Find total correlation energy
-      RS = ( 3 / (4*PI*DT) )**THD
-      KF = (3 * PI**2 * DT)**THD
-      KS = SQRT( 4 * KF / PI )
-      ZETA = ( D(1) - D(2) ) / DT
-      ZETA = MAX( -1.D0+DENMIN, ZETA )
-      ZETA = MIN(  1.D0-DENMIN, ZETA )
-      PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
-      T = GDMT / (2 * PHI * KS * DT)
-      F1 = ECUNIF / GAMMA / PHI**3
-      F2 = EXP(-F1)
-      A = BETA / GAMMA / (F2-1)
-      F3 = T**2 + A * T**4
-      F4 = BETA/GAMMA * F3 / (1 + A*F3)
-      H = GAMMA * PHI**3 * LOG( 1 + F4 )
-      FC = ECUNIF + H
-
-C Find correlation energy derivatives
-      DRSDD = - (THD * RS / DT)
-      DKFDD =   THD * KF / DT
-      DKSDD = HALF * KS * DKFDD / KF
-      DZDD(1) =   1 / DT - ZETA / DT
-      DZDD(2) = - (1 / DT) - ZETA / DT
-      DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
-      DO 40 IS = 1,2
-        DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
-        DPDD = DPDZ * DZDD(IS)
-        DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
-        DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
-        DF2DD = (- F2) * DF1DD
-        DADD = (- A) * DF2DD / (F2-1)
-        DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
-        DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
-        DHDD = 3 * H * DPDD / PHI
-        DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
-        DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
-
-        DO 30 IX = 1,3
-          DTDGD = (T / GDMT) * GDT(IX) / GDMT
-          DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
-          DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 
-          DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
-          DFCDGD(IX,IS) = DT * DHDGD
-   30   CONTINUE
-   40 CONTINUE
-
-C Find exchange energy and potential
-      FX = 0
-      DO 60 IS = 1,2
-        DS(IS)   = MAX( DENMIN, 2 * D(IS) )
-        GDMS = MAX( GDMIN, 2 * GDM(IS) )
-        KFS = (3 * PI**2 * DS(IS))**THD
-        S = GDMS / (2 * KFS * DS(IS))
-        F1 = 1 + MU * S**2 / KAPPA
-        F = 1 + KAPPA - KAPPA / F1
-c
-c       Note nspin=1 in call to exchng...
-c
-        CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
-        FX = FX + DS(IS) * EXUNIF * F
-
-        DKFDD = THD * KFS / DS(IS)
-        DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
-        DF1DD = 2 * (F1-1) * DSDD / S
-        DFDD = KAPPA * DF1DD / F1**2
-        DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
-
-        DO 50 IX = 1,3
-          GDS = 2 * GD(IX,IS)
-          DSDGD = (S / GDMS) * GDS / GDMS
-          DF1DGD = 2 * MU * S * DSDGD / KAPPA
-          DFDGD = KAPPA * DF1DGD / F1**2
-          DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
-   50   CONTINUE
-   60 CONTINUE
-      FX = HALF * FX / DT
-
-C Set output arguments
-      EX = FX
-      EC = FC
-      DO 90 IS = 1,nspin
-        DEXDD(IS) = DFXDD(IS)
-        DECDD(IS) = DFCDD(IS)
-        DO 80 IX = 1,3
-          DEXDGD(IX,IS) = DFXDGD(IX,IS)
-          DECDGD(IX,IS) = DFCDGD(IX,IS)
-   80   CONTINUE
-   90 CONTINUE
-
-      END SUBROUTINE PBEformXC
-
-
-
-      SUBROUTINE PBEXC( IREL, nspin, Dens, GDens,
-     .                  EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation.
-C Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
-C Modified to call PBEformXC by J.M.Soler. December 2009.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C EXCHNG, PW92C
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables
-      real(dp):: BETA, KAPPA, MU, PI
-
-C Fix values for PBE functional parameters
-      PI = 4 * ATAN(1._dp)
-      BETA = 0.066725_dp       ! From grad. exp. for correl. rs->0
-      MU = BETA * PI**2 / 3    ! From Jell. response for x+c
-      KAPPA = 0.804_dp         ! From general Lieb-Oxford bound
-
-C Call PBE routine with appropriate values for beta, mu, and kappa.
-      CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens,
-     .                EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-      END SUBROUTINE PBEXC
-
-
-
-      SUBROUTINE REVPBEXC( IREL, nspin, Dens, GDens,
-     .                     EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements revPBE: revised Perdew-Burke-Ernzerhof GGA.
-C Ref: Y. Zhang & W. Yang, Phys. Rev. Lett. 80, 890 (1998).
-C Same interface as PBEXC.
-C revPBE parameters introduced by E. Artacho in January 2006
-C Modified to call PBEformXC by J.M.Soler. December 2009.
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables
-      real(dp):: BETA, KAPPA, MU, PI
-
-C Fix values for PBE functional parameters
-      PI = 4 * ATAN(1._dp)
-      BETA = 0.066725_dp       ! From grad. exp. for correl. rs->0
-      MU = BETA * PI**2 / 3    ! From Jell. response for x+c
-      KAPPA = 1.245_dp         ! From fit of molecular energies
-
-C Call PBE routine with appropriate values for beta, mu, and kappa.
-      CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens,
-     .                EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-      END SUBROUTINE REVPBEXC
-
-
-
-      SUBROUTINE PBESOLXC( IREL, NSPIN, DENS, GDENS,
-     .                  EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation.
-C with the revised parameters for solids (PBEsol).
-C Ref: J.P.Perdew et al, PRL 100, 136406 (2008)
-C Same interface as PBEXC.
-C Modified by J.D. Gale for PBEsol. May 2009.
-C Modified to call PBEformXC by J.M.Soler. December 2009.
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables
-      real(dp):: BETA, KAPPA, MU, PI
-
-C Fix values for PBE functional parameters
-      PI = 4 * ATAN(1._dp)
-      BETA = 0.046_dp          ! From fit of Jell. surf. energies
-      MU = 10.0_dp / 81.0_dp   ! From grad. exp. for exchange.
-      KAPPA = 0.804_dp         ! From general Lieb-Oxford bound
-
-C Call PBE routine with appropriate values for beta, mu, and kappa.
-      CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens,
-     .                EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-      END SUBROUTINE PBESOLXC
-
-
-
-      SUBROUTINE PBEJsJrLOxc( IREL, NSPIN, DENS, GDENS,
-     .                      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
-C functional form, with revised parameters of Capelle et al:
-C   Js refers to Jellium surface energies, that fix parameter beta
-C   Jr refers to Jellium response, that fixes parameter mu
-C   LO refers to Lieb-Oxford bound, that fixes parameter kappa
-C Refs: L.S.Pedroza et al, PRB 79, 201106 (2009)
-C       M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009)
-C Same interface as PBEXC. J.M.Soler. December 2009.
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables
-      real(dp):: BETA, KAPPA, MU, PI
-
-C Fix values for PBE functional parameters
-      PI = 4 * ATAN(1._dp)
-      BETA = 0.046_dp          ! From fit of Jell. surf. energies
-      MU = BETA * PI**2 / 3    ! From Jell. response for x+c
-      KAPPA = 0.804_dp         ! From general Lieb-Oxford bound
-
-C Call PBE routine with appropriate values for beta, mu, and kappa.
-      CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens,
-     .                EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-      END SUBROUTINE PBEJsJrLOxc
-
-
-
-      SUBROUTINE PBEJsJrHEGxc( IREL, NSPIN, DENS, GDENS,
-     .                      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
-C functional form, with revised parameters of Capelle et al:
-C   Js refers to Jellium surface energies, that fix parameter beta
-C   Jr refers to Jellium response, that fixes parameter mu
-C   HGE refers to the Lieb-Oxford bound for the low-density limit of
-C       the homogeneous electron gas, that fixes parameter kappa
-C Refs: L.S.Pedroza et al, PRB 79, 201106 (2009)
-C       M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009)
-C Same interface as PBEXC. J.M.Soler. December 2009.
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables
-      real(dp):: BETA, KAPPA, MU, PI
-
-C Fix values for PBE functional parameters
-      PI = 4 * ATAN(1._dp)
-      BETA = 0.046_dp          ! From fit of Jell. surf. energies
-      MU = BETA * PI**2 / 3    ! From Jell. response for x+c
-      KAPPA = 0.552_dp         ! From Lieb-Oxford bound for HEG
-
-C Call PBE routine with appropriate values for beta, mu, and kappa.
-      CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens,
-     .                EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-      END SUBROUTINE PBEJsJrHEGxc
-
-
-
-      SUBROUTINE PBEGcGxLOxc( IREL, NSPIN, DENS, GDENS,
-     .                      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
-C functional form, with revised parameters of Capelle et al:
-C   Gc refers to gradient exp. for correl., that fixes parameter beta
-C   Gx refers to grad. expansion for exchange, that fixes parameter mu
-C   LO refers to Lieb-Oxford bound, that fixes parameter kappa
-C Refs: L.S.Pedroza et al, PRB 79, 201106 (2009)
-C       M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009)
-C Same interface as PBEXC. J.M.Soler. December 2009.
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables
-      real(dp):: BETA, KAPPA, MU, PI
-
-C Fix values for PBE functional parameters
-      PI = 4 * ATAN(1._dp)
-      BETA = 0.066725_dp       ! From grad. exp. for correl. rs->0
-      MU = 10._dp / 81._dp     ! From grad. exp. for exchange.
-      KAPPA = 0.804_dp         ! From general Lieb-Oxford bound
-
-C Call PBE routine with appropriate values for beta, mu, and kappa.
-      CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens,
-     .                EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-      END SUBROUTINE PBEGcGxLOxc
-
-
-
-      SUBROUTINE PBEGcGxHEGxc( IREL, NSPIN, DENS, GDENS,
-     .                      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
-C functional form, with revised parameters of Capelle et al:
-C   Gc refers to gradient exp. for correl., that fixes parameter beta
-C   Gx refers to grad. expansion for exchange, that fixes parameter mu
-C   HGE refers to the Lieb-Oxford bound for the low-density limit of
-C       the homogeneous electron gas, that fixes parameter kappa
-C Refs: L.S.Pedroza et al, PRB 79, 201106 (2009)
-C       M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009)
-C Same interface as PBEXC. J.M.Soler. December 2009.
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables
-      real(dp):: BETA, KAPPA, MU, PI
-
-C Fix values for PBE functional parameters
-      PI = 4 * ATAN(1._dp)
-      BETA = 0.066725_dp       ! From grad. exp. for correl. rs->0
-      MU = 10._dp / 81._dp     ! From grad. exp. for exchange.
-      KAPPA = 0.552_dp         ! From Lieb-Oxford bound for HEG
-
-C Call PBE routine with appropriate values for beta, mu, and kappa.
-      CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens,
-     .                EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-      END SUBROUTINE PBEGcGxHEGxc
-
-
-
-      SUBROUTINE PW91XC( IREL, nspin, Dens, GDens,
-     .                  EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Wang91 Generalized-Gradient-Approximation.
-C Ref: JCP 100, 1290 (1994)
-C Written by J.L. Martins  August 2000
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C EXCHNG, PW92C
-C ********************************************************************
-
-      implicit          none
-      INTEGER           IREL, nspin
-      real(dp)          Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
-     .                  DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
-
-C Internal variables
-      INTEGER
-     .  IS, IX
-      real(dp)
-     .  A, BETA, D(2), DADD, DECUDD, DENMIN, 
-     .  DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
-     .  DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
-     .  DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, 
-     .  DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), 
-     .  EC, ECUNIF, EX, EXUNIF,
-     .  F, F1, F2, F3, F4, FC, FX, FOUTHD,
-     .  GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
-     .  H, HALF, KF, KFS, KS, PHI, PI, RS, S,
-     .  T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
-     
-      real(dp)          F5, F6, F7, F8, ASINHS
-      real(dp)          DF5DD,DF6DD,DF7DD,DF8DD
-      real(dp)          DF1DS, DF2DS, DF3DS, DFDS, DF7DGD
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      PARAMETER ( DENMIN = 1.D-12 )
-      PARAMETER ( GDMIN  = 1.D-12 )
-
-C Fix some numerical parameters
-      PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
-     .            THD=1.D0/3.D0, THRHLF=1.5D0,
-     .            TWO=2.D0, TWOTHD=2.D0/3.D0 )
-
-C Fix some more numerical constants
-      PI = 4.0_dp * ATAN(1.0_dp)
-      BETA = 15.75592_dp * 0.004235_dp
-      GAMMA = BETA**2 / (2.0_dp * 0.09_dp)
-
-C Translate density and its gradient to new variables
-      IF (nspin .EQ. 1) THEN
-        D(1) = HALF * Dens(1)
-        D(2) = D(1)
-        DT = MAX( DENMIN, Dens(1) )
-        DO 10 IX = 1,3
-          GD(IX,1) = HALF * GDens(IX,1)
-          GD(IX,2) = GD(IX,1)
-          GDT(IX) = GDens(IX,1)
-   10   CONTINUE
-      ELSE
-        D(1) = Dens(1)
-        D(2) = Dens(2)
-        DT = MAX( DENMIN, Dens(1)+Dens(2) )
-        DO 20 IX = 1,3
-          GD(IX,1) = GDens(IX,1)
-          GD(IX,2) = GDens(IX,2)
-          GDT(IX) = GDens(IX,1) + GDens(IX,2)
-   20   CONTINUE
-      ENDIF
-      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = MAX( GDMIN, GDMT )
-
-C Find local correlation energy and potential
-      CALL PW92C( 2, D, ECUNIF, VCUNIF )
-
-C Find total correlation energy
-      RS = ( 3 / (4*PI*DT) )**THD
-      KF = (3 * PI**2 * DT)**THD
-      KS = SQRT( 4 * KF / PI )
-      S = GDMT / (2 * KF * DT)
-      T = GDMT / (2 * KS * DT)
-      ZETA = ( D(1) - D(2) ) / DT
-      ZETA = MAX( -1.D0+DENMIN, ZETA )
-      ZETA = MIN(  1.D0-DENMIN, ZETA )
-      PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
-      F1 = ECUNIF / GAMMA / PHI**3
-      F2 = EXP(-F1)
-      A = BETA / GAMMA / (F2-1)
-      F3 = T**2 + A * T**4
-      F4 = BETA/GAMMA * F3 / (1 + A*F3)
-      F5 = 0.002568D0 + 0.023266D0*RS + 7.389D-6*RS**2
-      F6 = 1.0D0 + 8.723D0*RS + 0.472D0*RS**2 + 0.07389D0*RS**3
-      F7 = EXP(-100.0D0 * S**2 * PHI**4)
-      F8 =  15.75592D0*(0.001667212D0 + F5/F6 -0.004235D0 + 
-     .          3.0D0*0.001667212D0/7.0D0)
-      H = GAMMA * PHI**3 * LOG( 1 + F4 ) + F8 * T**2 * F7
-      FC = ECUNIF + H
-
-C Find correlation energy derivatives
-      DRSDD = - THD * RS / DT
-      DKFDD =   THD * KF / DT
-      DKSDD = HALF * KS * DKFDD / KF
-      DZDD(1) =   1 / DT - ZETA / DT
-      DZDD(2) = - 1 / DT - ZETA / DT
-      DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
-      DO 40 IS = 1,2
-        DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
-        DPDD = DPDZ * DZDD(IS)
-        DSDD = -S*DKFDD/KF - S/DT    ! JMS: corrected, May.2014
-        DTDD = -T*DKSDD/KS - T/DT
-        DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
-        DF2DD = - F2 * DF1DD
-        DADD = - A * DF2DD / (F2-1)
-        DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
-        DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
-        DF5DD = (0.023266D0 + 2.0D0*7.389D-6*RS)*DRSDD
-        DF6DD = (8.723D0 + 2.0D0*0.472D0*RS
-     .            + 3.0D0*0.07389D0*RS**2)*DRSDD
-        DF7DD = -200.0D0 * S * PHI**4 * DSDD * F7
-     .         -100.0D0 * S**2 * 4.0D0* PHI**3 * DPDD * F7
-        DF8DD = 15.75592D0 * DF5DD/F6 - 15.75592D0*F5*DF6DD / F6**2
-        DHDD = 3 * H * DPDD / PHI
-        DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
-        DHDD = DHDD + DF8DD * T**2 * F7
-        DHDD = DHDD + F8 * 2*T*DTDD *F7
-        DHDD = DHDD + F8 * T**2 * DF7DD
-        
-        DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
-        DO 30 IX = 1,3
-          DTDGD = (T / GDMT) * GDT(IX) / GDMT
-          DSDGD = (S / GDMT) * GDT(IX) / GDMT
-          DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
-          DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 
-          DF7DGD = -200.0D0 * S * PHI**4 * DSDGD * F7
-          DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
-          DHDGD = DHDGD + F8 * 2*T*DTDGD *F7 + F8 * T**2 *DF7DGD
-          DFCDGD(IX,IS) = DT * DHDGD
-   30   CONTINUE
-   40 CONTINUE
-
-C Find exchange energy and potential
-      FX = 0
-      DO 60 IS = 1,2
-        DS(1) = MAX( DENMIN, 2 * D(IS) )
-        GDMS = MAX( GDMIN, 2 * GDM(IS) )
-        KFS = (3 * PI**2 * DS(1))**THD
-        S = GDMS / (2 * KFS * DS(1))
-        F4 = SQRT(1.0D0 + (7.7956D0*S)**2)
-        ASINHS = LOG(7.7956D0*S + F4)
-        F1 = 1.0D0 + 0.19645D0 * S * ASINHS
-        F2 = 0.2743D0 - 0.15084D0*EXP(-100.0D0*S*S)
-        F3 = 1.0D0 / (F1 + 0.004D0 * S*S*S*S)
-        F = (F1 + F2 * S*S ) * F3
-     .       
-        CALL EXCHNG( IREL, 1, DS, EXUNIF, VXUNIF )
-        FX = FX + DS(1) * EXUNIF * F
-
-        DKFDD = THD * KFS / DS(1)
-        DSDD = S * ( -DKFDD/KFS - 1/DS(1) )
-        DF1DS = 0.19645D0 * ASINHS +
-     .    0.19645D0 * S * 7.7956D0 / F4
-        DF2DS = 0.15084D0*200.0D0*S*EXP(-100.0D0*S*S)
-        DF3DS = - F3*F3 * (DF1DS + 4.0D0*0.004D0 * S*S*S)
-        DFDS =  DF1DS * F3 + DF2DS * S*S * F3 + 2.0D0 * S * F2 * F3
-     .            + (F1 + F2 * S*S ) * DF3DS   
-        DFXDD(IS) = VXUNIF(1) * F + DS(1) * EXUNIF * DFDS * DSDD
-
-        DO 50 IX = 1,3
-          GDS = 2 * GD(IX,IS)
-          DSDGD = (S / GDMS) * GDS / GDMS
-          DFDGD = DFDS * DSDGD
-          DFXDGD(IX,IS) = DS(1) * EXUNIF * DFDGD
-   50   CONTINUE
-   60 CONTINUE
-      FX = HALF * FX / DT
-
-C Set output arguments
-      EX = FX
-      EC = FC
-      DO 90 IS = 1,nspin
-        DEXDD(IS) = DFXDD(IS)
-        DECDD(IS) = DFCDD(IS)
-        DO 80 IX = 1,3
-          DEXDGD(IX,IS) = DFXDGD(IX,IS)
-          DECDGD(IX,IS) = DFCDGD(IX,IS)
-   80   CONTINUE
-   90 CONTINUE
-
-      END SUBROUTINE PW91XC
-
-
-
-       subroutine blypxc(nspin,dens,gdens,EX,EC,
-     .                   dEXdd,dECdd,dEXdgd,dECdgd) 
-c ***************************************************************
-c Implements Becke gradient exchange functional (A.D. 
-c Becke, Phys. Rev. A 38, 3098 (1988)) and Lee, Yang, Parr
-c correlation functional (C. Lee, W. Yang, R.G. Parr, Phys. Rev. B
-c 37, 785 (1988)), as modificated by Miehlich,Savin,Stoll and Preuss,
-c Chem. Phys. Lett. 157,200 (1989). See also Johnson, Gill and Pople,
-c J. Chem. Phys. 98, 5612 (1993). Some errors were detected in this
-c last paper, so not all of the expressions correspond exactly to those
-c implemented here.
-c Written by Maider Machado. July 1998.
-c **************** INPUT ******************************************** 
-c integer nspin          : Number of spin polarizations (1 or 2)
-c real*8  dens(nspin)    : Total electron density (if nspin=1) or
-c                           spin electron density (if nspin=2)
-c real*8  gdens(3,nspin) : Total or spin density gradient
-c ******** OUTPUT *****************************************************
-c real*8  ex             : Exchange energy density
-c real*8  ec             : Correlation energy density
-c real*8  dexdd(nspin)   : Partial derivative
-c                           d(DensTot*Ex)/dDens(ispin),
-c                           where DensTot = Sum_ispin( Dens(ispin) )
-c                          For a constant density, this is the
-c                          exchange potential
-c real*8  decdd(nspin)   : Partial derivative
-c                           d(DensTot*Ec)/dDens(ispin),
-c                           where DensTot = Sum_ispin( Dens(ispin) )
-c                          For a constant density, this is the
-c                          correlation potential
-c real*8  dexdgd(3,nspin): Partial derivative
-c                           d(DensTot*Ex)/d(GradDens(i,ispin))
-c real*8  decdgd(3,nspin): Partial derivative
-c                           d(DensTot*Ec)/d(GradDens(i,ispin))
-c ********* UNITS ****************************************************
-c Lengths in Bohr
-c Densities in electrons per Bohr**3
-c Energies in Hartrees
-c Gradient vectors in cartesian coordinates
-c ********************************************************************
- 
-      implicit none
-
-      integer nspin
-      real(dp)   dens(nspin), gdens(3,nspin), EX, EC,
-     .           dEXdd(nspin), dECdd(nspin), dEXdgd(3,nspin),
-     .           dECdgd(3,nspin)
-
-c Internal variables
-      integer is,ix
-      real(dp)   pi, beta, thd, tthd, thrhlf, half, fothd,
-     .           d(2),gd(3,2),dmin, ash,gdm(2),denmin,dt, 
-     .           g(2),x(2),a,b,c,dd,onzthd,gdmin,     
-     .           ga, gb, gc,becke,dbecgd(3,2),
-     .           dgdx(2), dgdxa, dgdxb, dgdxc,dgdxd,dbecdd(2),
-     .           den,omega, domega, delta, ddelta,cf,
-     .           gam11, gam12, gam22, LYPa, LYPb1,
-     .           LYPb2,dLYP11,dLYP12,dLYP22,LYP,
-     .           dd1g11,dd1g12,dd1g22,dd2g12,dd2g11,dd2g22,
-     .           dLYPdd(2),dg11dd(3,2),dg22dd(3,2),
-     .           dg12dd(3,2),dLYPgd(3,2)
-  
-c Lower bounds of density and its gradient to avoid divisions by zero
-      parameter ( denmin=1.0e-8_dp )
-      parameter ( gdmin=1.0e-8_dp )
-      parameter ( dmin=1.0e-5_dp ) 
-
-c Fix some numerical parameters 
-      parameter ( thd = 1.d0/3.d0, tthd=2.d0/3.d0 )
-      parameter ( thrhlf=1.5d0, half=0.5d0,
-     .            fothd=4.d0/3.d0, onzthd=11.d0/3.d0)
-
-c Empirical parameter for Becke exchange functional (a.u.)
-      parameter(beta= 0.0042d0) 
-
-c Constants for LYP functional (a.u.) 
-      parameter(a=0.04918d0, b=0.132d0, c=0.2533d0, dd=0.349d0)
-
-      pi= 4*atan(1.d0)
-
-c Translate density and its gradient to new variables
-      if (nspin .eq. 1) then
-        d(1) = half * dens(1)
-        d(1) = max(denmin,d(1))
-        d(2) = d(1)
-        dt = max( denmin, dens(1) )
-        do ix = 1,3
-          gd(ix,1) = half * gdens(ix,1)    
-          gd(ix,2) = gd(ix,1)
-        enddo 
-      else
-        d(1) = dens(1)
-        d(2) = dens(2)
-        do is=1,2
-         d(is) = max (denmin,d(is))
-        enddo
-        dt = max( denmin, dens(1)+dens(2) )  
-        do ix = 1,3
-          gd(ix,1) = gdens(ix,1)
-          gd(ix,2) = gdens(ix,2)
-        enddo
       endif
-
-      gdm(1) = sqrt( gd(1,1)**2 + gd(2,1)**2 + gd(3,1)**2 )
-      gdm(2) = sqrt( gd(1,2)**2 + gd(2,2)**2 + gd(3,2)**2 )
- 
+    else   ! Collinear spin => just copy derivatives to output arrays
+      dEXdD(1:nSpin) = dEXdDD(1:nSpin)
+      dECdD(1:nSpin) = dECdDD(1:nSpin)
+      dEXdGD(:,1:nSpin) = dEXdGDD(:,1:nSpin)
+      dECdGD(:,1:nSpin) = dECdGDD(:,1:nSpin)
+    end if ! (nSpin==4)
+
+  END SUBROUTINE GGAXC
+
+  
+  
+  SUBROUTINE PBEformXC( beta, mu, kappa, iRel, nSpin, Dens, GDens, &
+      EX, EC, dEXdD, dECdD, dEXdGD, dECdGD )
+
+    ! *********************************************************************
+    ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
+    ! functional form, but with variable values for parameters beta, mu, and
+    ! kappa. Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
+    ! Written by L.C.Balbas and J.M.Soler. December 1996. 
+    ! ******** INPUT ******************************************************
+    ! REAL*8  BETA           : Parameter beta of the PBE functional
+    ! REAL*8  MU             : Parameter mu of the PBE functional
+    ! REAL*8  KAPPA          : Parameter kappa of the PBE functional
+    ! INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
+    ! INTEGER nspin          : Number of spin polarizations (1 or 2)
+    ! REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
+    !                           spin electron density (if nspin=2)
+    ! REAL*8  GDens(3,nspin) : Total or spin density gradient
+    ! ******** OUTPUT *****************************************************
+    ! REAL*8  EX             : Exchange energy density
+    ! REAL*8  EC             : Correlation energy density
+    ! REAL*8  DEXDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ex)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          exchange potential
+    ! REAL*8  DECDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ec)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          correlation potential
+    ! REAL*8  DEXDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ex)/d(GradDens(i,ispin))
+    ! REAL*8  DECDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ec)/d(GradDens(i,ispin))
+    ! ********* UNITS ****************************************************
+    ! Lengths in Bohr
+    ! Densities in electrons per Bohr**3
+    ! Energies in Hartrees
+    ! Gradient vectors in cartesian coordinates
+    ! ********* ROUTINES CALLED ******************************************
+    ! EXCHNG, PW92C
+    ! ********************************************************************
+
+    ! Input
+    real(dp),intent(in) :: beta           ! Parameter of PBE functional
+    real(dp),intent(in) :: mu             ! Parameter of PBE functional
+    real(dp),intent(in) :: kappa          ! Parameter of PBE functional
+    integer, intent(in) :: iRel           ! Relativistic exchange? 0=no, 1=yes
+    integer, intent(in) :: nSpin          ! Number of spin components
+    real(dp),intent(in) :: Dens(nSpin)    ! Local electron (spin) density
+    real(dp),intent(in) :: GDens(3,nSpin) ! Gradient of electron density
+
+    ! Output
+    real(dp),intent(out):: EX             ! Exchange energy per electron
+    real(dp),intent(out):: EC             ! Correlation energy per electron
+    real(dp),intent(out):: dEXdD(nSpin)   ! dEx/dDens, Ex=Int(dens*epsX)
+    real(dp),intent(out):: dECdD(nSpin)   ! dEc/dDens
+    real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens)
+    real(dp),intent(out):: dECdGD(3,nSpin) ! dEc/dGrad(Dens)
+
+    ! Internal variables
+    INTEGER :: IS, IX
+    real(dp) &
+        A, D(2), DADD, DECUDD, DENMIN, &
+        DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD, &
+        DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2), &
+        DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, &
+        DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), &
+        ECUNIF, EXUNIF, &
+        F, F1, F2, F3, F4, FC, FX, FOUTHD, &
+        GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), &
+        H, HALF, KF, KFS, KS, PHI, PI, RS, S, &
+        T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
+
+    ! Lower bounds of density and its gradient to avoid divisions by zero
+    PARAMETER ( DENMIN = 1.D-12 )
+    PARAMETER ( GDMIN  = 1.D-12 )
+
+    ! Fix some numerical parameters
+    PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, &
+        THD=1.D0/3.D0, THRHLF=1.5D0, &
+        TWO=2.D0, TWOTHD=2.D0/3.D0 )
+
+    ! Fix some more numerical constants
+    PI = 4 * ATAN(1.D0)
+    GAMMA = (1 - LOG(TWO)) / PI**2
+
+    ! Translate density and its gradient to new variables
+    IF (nspin .EQ. 1) THEN
+      D(1) = HALF * Dens(1)
+      D(2) = D(1)
+      DT = MAX( DENMIN, Dens(1) )
+      DO IX = 1,3
+        GD(IX,1) = HALF * GDens(IX,1)
+        GD(IX,2) = GD(IX,1)
+        GDT(IX) = GDens(IX,1)
+      end DO
+    ELSE
+      D(1) = Dens(1)
+      D(2) = Dens(2)
+      DT = MAX( DENMIN, Dens(1)+Dens(2) )
+      DO IX = 1,3
+        GD(IX,1) = GDens(IX,1)
+        GD(IX,2) = GDens(IX,2)
+        GDT(IX) = GDens(IX,1) + GDens(IX,2)
+      end DO
+    ENDIF
+    GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
+    GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
+    GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
+    GDMT = MAX( GDMIN, GDMT )
+
+    ! Find local correlation energy and potential
+    CALL PW92C( 2, D, ECUNIF, VCUNIF )
+
+    ! Find total correlation energy
+    RS = ( 3 / (4*PI*DT) )**THD
+    KF = (3 * PI**2 * DT)**THD
+    KS = SQRT( 4 * KF / PI )
+    ZETA = ( D(1) - D(2) ) / DT
+    ZETA = MAX( -1.D0+DENMIN, ZETA )
+    ZETA = MIN(  1.D0-DENMIN, ZETA )
+    PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
+    T = GDMT / (2 * PHI * KS * DT)
+    F1 = ECUNIF / GAMMA / PHI**3
+    F2 = EXP(-F1)
+    A = BETA / GAMMA / (F2-1)
+    F3 = T**2 + A * T**4
+    F4 = BETA/GAMMA * F3 / (1 + A*F3)
+    H = GAMMA * PHI**3 * LOG( 1 + F4 )
+    FC = ECUNIF + H
+
+    ! Find correlation energy derivatives
+    DRSDD = - (THD * RS / DT)
+    DKFDD =   THD * KF / DT
+    DKSDD = HALF * KS * DKFDD / KF
+    DZDD(1) =   1 / DT - ZETA / DT
+    DZDD(2) = - (1 / DT) - ZETA / DT
+    DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
+    DO IS = 1,2
+      DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
+      DPDD = DPDZ * DZDD(IS)
+      DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
+      DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
+      DF2DD = (- F2) * DF1DD
+      DADD = (- A) * DF2DD / (F2-1)
+      DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
+      DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
+      DHDD = 3 * H * DPDD / PHI
+      DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
+      DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
+
+      DO IX = 1,3
+        DTDGD = (T / GDMT) * GDT(IX) / GDMT
+        DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
+        DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 
+        DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
+        DFCDGD(IX,IS) = DT * DHDGD
+      end DO
+    end DO
+
+    ! Find exchange energy and potential
+    FX = 0
+    DO IS = 1,2
+      DS(IS)   = MAX( DENMIN, 2 * D(IS) )
+      GDMS = MAX( GDMIN, 2 * GDM(IS) )
+      KFS = (3 * PI**2 * DS(IS))**THD
+      S = GDMS / (2 * KFS * DS(IS))
+      F1 = 1 + MU * S**2 / KAPPA
+      F = 1 + KAPPA - KAPPA / F1
+      !
+      !       Note nspin=1 in call to exchng...
+      !
+      CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
+      FX = FX + DS(IS) * EXUNIF * F
+
+      DKFDD = THD * KFS / DS(IS)
+      DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
+      DF1DD = 2 * (F1-1) * DSDD / S
+      DFDD = KAPPA * DF1DD / F1**2
+      DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
+
+      DO IX = 1,3
+        GDS = 2 * GD(IX,IS)
+        DSDGD = (S / GDMS) * GDS / GDMS
+        DF1DGD = 2 * MU * S * DSDGD / KAPPA
+        DFDGD = KAPPA * DF1DGD / F1**2
+        DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
+      end DO
+    end DO
+    FX = HALF * FX / DT
+
+    ! Set output arguments
+    EX = FX
+    EC = FC
+    DO IS = 1,nspin
+      DEXDD(IS) = DFXDD(IS)
+      DECDD(IS) = DFCDD(IS)
+      DO IX = 1,3
+        DEXDGD(IX,IS) = DFXDGD(IX,IS)
+        DECDGD(IX,IS) = DFCDGD(IX,IS)
+      end DO
+    end DO
+  END SUBROUTINE PBEformXC
+
+
+
+  SUBROUTINE PBEXC( IREL, nspin, Dens, GDens, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation.
+    ! Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
+    ! Modified to call PBEformXC by J.M.Soler. December 2009.
+    ! ******** INPUT ******************************************************
+    ! INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
+    ! INTEGER nspin          : Number of spin polarizations (1 or 2)
+    ! REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
+    !                           spin electron density (if nspin=2)
+    ! REAL*8  GDens(3,nspin) : Total or spin density gradient
+    ! ******** OUTPUT *****************************************************
+    ! REAL*8  EX             : Exchange energy density
+    ! REAL*8  EC             : Correlation energy density
+    ! REAL*8  DEXDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ex)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          exchange potential
+    ! REAL*8  DECDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ec)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          correlation potential
+    ! REAL*8  DEXDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ex)/d(GradDens(i,ispin))
+    ! REAL*8  DECDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ec)/d(GradDens(i,ispin))
+    ! ********* UNITS ****************************************************
+    ! Lengths in Bohr
+    ! Densities in electrons per Bohr**3
+    ! Energies in Hartrees
+    ! Gradient vectors in cartesian coordinates
+    ! ********* ROUTINES CALLED ******************************************
+    ! EXCHNG, PW92C
+    ! ********************************************************************
+
+    ! Passed arguments
+    integer, intent(in) :: IREL, NSPIN
+    real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
+    real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), &
+        DEXDD(NSPIN), DEXDGD(3,NSPIN)
+
+    ! Internal variables
+    real(dp):: BETA, KAPPA, MU, PI
+
+    ! Fix values for PBE functional parameters
+    PI = 4 * ATAN(1._dp)
+    BETA = 0.066725_dp       ! From grad. exp. for correl. rs->0
+    MU = BETA * PI**2 / 3    ! From Jell. response for x+c
+    KAPPA = 0.804_dp         ! From general Lieb-Oxford bound
+
+    ! Call PBE routine with appropriate values for beta, mu, and kappa.
+    CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, &
+        EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+  END SUBROUTINE PBEXC
+
+
+
+  SUBROUTINE REVPBEXC( IREL, nspin, Dens, GDens,&
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements revPBE: revised Perdew-Burke-Ernzerhof GGA.
+    ! Ref: Y. Zhang & W. Yang, Phys. Rev. Lett. 80, 890 (1998).
+    ! Same interface as PBEXC.
+    ! revPBE parameters introduced by E. Artacho in January 2006
+    ! Modified to call PBEformXC by J.M.Soler. December 2009.
+    ! ********************************************************************
+
+    ! Passed arguments
+    integer, intent(in) :: IREL, NSPIN
+    real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
+    real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), &
+        DEXDD(NSPIN), DEXDGD(3,NSPIN)
+
+    ! Internal variables
+    real(dp):: BETA, KAPPA, MU, PI
+
+    ! Fix values for PBE functional parameters
+    PI = 4 * ATAN(1._dp)
+    BETA = 0.066725_dp       ! From grad. exp. for correl. rs->0
+    MU = BETA * PI**2 / 3    ! From Jell. response for x+c
+    KAPPA = 1.245_dp         ! From fit of molecular energies
+
+    ! Call PBE routine with appropriate values for beta, mu, and kappa.
+    CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, &
+        EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+  END SUBROUTINE REVPBEXC
+
+
+
+  SUBROUTINE PBESOLXC( IREL, NSPIN, DENS, GDENS, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation.
+    ! with the revised parameters for solids (PBEsol).
+    ! Ref: J.P.Perdew et al, PRL 100, 136406 (2008)
+    ! Same interface as PBEXC.
+    ! Modified by J.D. Gale for PBEsol. May 2009.
+    ! Modified to call PBEformXC by J.M.Soler. December 2009.
+    ! ********************************************************************
+
+    ! Passed arguments
+    integer, intent(in) :: IREL, NSPIN
+    real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
+    real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), &
+        DEXDD(NSPIN), DEXDGD(3,NSPIN)
+
+    ! Internal variables
+    real(dp):: BETA, KAPPA, MU, PI
+
+    ! Fix values for PBE functional parameters
+    PI = 4 * ATAN(1._dp)
+    BETA = 0.046_dp          ! From fit of Jell. surf. energies
+    MU = 10.0_dp / 81.0_dp   ! From grad. exp. for exchange.
+    KAPPA = 0.804_dp         ! From general Lieb-Oxford bound
+
+    ! Call PBE routine with appropriate values for beta, mu, and kappa.
+    CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, &
+        EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+  END SUBROUTINE PBESOLXC
+
+
+
+  SUBROUTINE PBEJsJrLOxc( IREL, NSPIN, DENS, GDENS, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
+    ! functional form, with revised parameters of Capelle et al:
+    !   Js refers to Jellium surface energies, that fix parameter beta
+    !   Jr refers to Jellium response, that fixes parameter mu
+    !   LO refers to Lieb-Oxford bound, that fixes parameter kappa
+    ! Refs: L.S.Pedroza et al, PRB 79, 201106 (2009)
+    !       M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009)
+    ! Same interface as PBEXC. J.M.Soler. December 2009.
+    ! ********************************************************************
+
+    ! Passed arguments
+    integer, intent(in) :: IREL, NSPIN
+    real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
+    real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), &
+        DEXDD(NSPIN), DEXDGD(3,NSPIN)
+
+    ! Internal variables
+    real(dp):: BETA, KAPPA, MU, PI
+
+    ! Fix values for PBE functional parameters
+    PI = 4 * ATAN(1._dp)
+    BETA = 0.046_dp          ! From fit of Jell. surf. energies
+    MU = BETA * PI**2 / 3    ! From Jell. response for x+c
+    KAPPA = 0.804_dp         ! From general Lieb-Oxford bound
+
+    ! Call PBE routine with appropriate values for beta, mu, and kappa.
+    CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, &
+        EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+  END SUBROUTINE PBEJsJrLOxc
+
+
+
+  SUBROUTINE PBEJsJrHEGxc( IREL, NSPIN, DENS, GDENS, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
+    ! functional form, with revised parameters of Capelle et al:
+    !   Js refers to Jellium surface energies, that fix parameter beta
+    !   Jr refers to Jellium response, that fixes parameter mu
+    !   HGE refers to the Lieb-Oxford bound for the low-density limit of
+    !       the homogeneous electron gas, that fixes parameter kappa
+    ! Refs: L.S.Pedroza et al, PRB 79, 201106 (2009)
+    !       M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009)
+    ! Same interface as PBEXC. J.M.Soler. December 2009.
+    ! ********************************************************************
+
+    ! Passed arguments
+    integer, intent(in) :: IREL, NSPIN
+    real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
+    real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), &
+        DEXDD(NSPIN), DEXDGD(3,NSPIN)
+
+    ! Internal variables
+    real(dp):: BETA, KAPPA, MU, PI
+
+    ! Fix values for PBE functional parameters
+    PI = 4 * ATAN(1._dp)
+    BETA = 0.046_dp          ! From fit of Jell. surf. energies
+    MU = BETA * PI**2 / 3    ! From Jell. response for x+c
+    KAPPA = 0.552_dp         ! From Lieb-Oxford bound for HEG
+
+    ! Call PBE routine with appropriate values for beta, mu, and kappa.
+    CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, &
+        EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+  END SUBROUTINE PBEJsJrHEGxc
+
+
+
+  SUBROUTINE PBEGcGxLOxc( IREL, NSPIN, DENS, GDENS, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
+    ! functional form, with revised parameters of Capelle et al:
+    !   Gc refers to gradient exp. for correl., that fixes parameter beta
+    !   Gx refers to grad. expansion for exchange, that fixes parameter mu
+    !   LO refers to Lieb-Oxford bound, that fixes parameter kappa
+    ! Refs: L.S.Pedroza et al, PRB 79, 201106 (2009)
+    !       M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009)
+    ! Same interface as PBEXC. J.M.Soler. December 2009.
+    ! ********************************************************************
+
+    ! Passed arguments
+    integer, intent(in) :: IREL, NSPIN
+    real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
+    real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), &
+        DEXDD(NSPIN), DEXDGD(3,NSPIN)
+
+    ! Internal variables
+    real(dp):: BETA, KAPPA, MU, PI
+
+    ! Fix values for PBE functional parameters
+    PI = 4 * ATAN(1._dp)
+    BETA = 0.066725_dp       ! From grad. exp. for correl. rs->0
+    MU = 10._dp / 81._dp     ! From grad. exp. for exchange.
+    KAPPA = 0.804_dp         ! From general Lieb-Oxford bound
+
+    ! Call PBE routine with appropriate values for beta, mu, and kappa.
+    CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, &
+        EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+  END SUBROUTINE PBEGcGxLOxc
+
+
+
+  SUBROUTINE PBEGcGxHEGxc( IREL, NSPIN, DENS, GDENS, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation
+    ! functional form, with revised parameters of Capelle et al:
+    !   Gc refers to gradient exp. for correl., that fixes parameter beta
+    !   Gx refers to grad. expansion for exchange, that fixes parameter mu
+    !   HGE refers to the Lieb-Oxford bound for the low-density limit of
+    !       the homogeneous electron gas, that fixes parameter kappa
+    ! Refs: L.S.Pedroza et al, PRB 79, 201106 (2009)
+    !       M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009)
+    ! Same interface as PBEXC. J.M.Soler. December 2009.
+    ! ********************************************************************
+
+    ! Passed arguments
+    integer, intent(in) :: IREL, NSPIN
+    real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
+    real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), &
+        DEXDD(NSPIN), DEXDGD(3,NSPIN)
+
+    ! Internal variables
+    real(dp):: BETA, KAPPA, MU, PI
+
+    ! Fix values for PBE functional parameters
+    PI = 4 * ATAN(1._dp)
+    BETA = 0.066725_dp       ! From grad. exp. for correl. rs->0
+    MU = 10._dp / 81._dp     ! From grad. exp. for exchange.
+    KAPPA = 0.552_dp         ! From Lieb-Oxford bound for HEG
+
+    ! Call PBE routine with appropriate values for beta, mu, and kappa.
+    CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, &
+        EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+  END SUBROUTINE PBEGcGxHEGxc
+
+
+
+  SUBROUTINE PW91XC( IREL, nspin, Dens, GDens, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Perdew-Wang91 Generalized-Gradient-Approximation.
+    ! Ref: JCP 100, 1290 (1994)
+    ! Written by J.L. Martins  August 2000
+    ! ******** INPUT ******************************************************
+    ! INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
+    ! INTEGER nspin          : Number of spin polarizations (1 or 2)
+    ! REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
+    !                           spin electron density (if nspin=2)
+    ! REAL*8  GDens(3,nspin) : Total or spin density gradient
+    ! ******** OUTPUT *****************************************************
+    ! REAL*8  EX             : Exchange energy density
+    ! REAL*8  EC             : Correlation energy density
+    ! REAL*8  DEXDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ex)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          exchange potential
+    ! REAL*8  DECDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ec)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          correlation potential
+    ! REAL*8  DEXDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ex)/d(GradDens(i,ispin))
+    ! REAL*8  DECDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ec)/d(GradDens(i,ispin))
+    ! ********* UNITS ****************************************************
+    ! Lengths in Bohr
+    ! Densities in electrons per Bohr**3
+    ! Energies in Hartrees
+    ! Gradient vectors in cartesian coordinates
+    ! ********* ROUTINES CALLED ******************************************
+    ! EXCHNG, PW92C
+    ! ********************************************************************
+
+    INTEGER           IREL, nspin
+    real(dp)          Dens(nspin), DECDD(nspin), DECDGD(3,nspin), &
+        DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
+
+    ! Internal variables
+    INTEGER :: IS, IX
+    real(dp) &
+        A, BETA, D(2), DADD, DECUDD, DENMIN, &
+        DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD, &
+        DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2), &
+        DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, &
+        DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), &
+        EC, ECUNIF, EX, EXUNIF, &
+        F, F1, F2, F3, F4, FC, FX, FOUTHD, &
+        GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), &
+        H, HALF, KF, KFS, KS, PHI, PI, RS, S, &
+        T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
+
+    real(dp)          F5, F6, F7, F8, ASINHS
+    real(dp)          DF5DD,DF6DD,DF7DD,DF8DD
+    real(dp)          DF1DS, DF2DS, DF3DS, DFDS, DF7DGD
+
+    ! Lower bounds of density and its gradient to avoid divisions by zero
+    PARAMETER ( DENMIN = 1.D-12 )
+    PARAMETER ( GDMIN  = 1.D-12 )
+
+    ! Fix some numerical parameters
+    PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, &
+        THD=1.D0/3.D0, THRHLF=1.5D0, &
+        TWO=2.D0, TWOTHD=2.D0/3.D0 )
+
+    ! Fix some more numerical constants
+    PI = 4.0_dp * ATAN(1.0_dp)
+    BETA = 15.75592_dp * 0.004235_dp
+    GAMMA = BETA**2 / (2.0_dp * 0.09_dp)
+
+    ! Translate density and its gradient to new variables
+    IF (nspin .EQ. 1) THEN
+      D(1) = HALF * Dens(1)
+      D(2) = D(1)
+      DT = MAX( DENMIN, Dens(1) )
+      DO IX = 1,3
+        GD(IX,1) = HALF * GDens(IX,1)
+        GD(IX,2) = GD(IX,1)
+        GDT(IX) = GDens(IX,1)
+      end DO
+    ELSE
+      D(1) = Dens(1)
+      D(2) = Dens(2)
+      DT = MAX( DENMIN, Dens(1)+Dens(2) )
+      DO IX = 1,3
+        GD(IX,1) = GDens(IX,1)
+        GD(IX,2) = GDens(IX,2)
+        GDT(IX) = GDens(IX,1) + GDens(IX,2)
+      end DO
+    ENDIF
+    GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
+    GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
+    GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
+    GDMT = MAX( GDMIN, GDMT )
+
+    ! Find local correlation energy and potential
+    CALL PW92C( 2, D, ECUNIF, VCUNIF )
+
+    ! Find total correlation energy
+    RS = ( 3 / (4*PI*DT) )**THD
+    KF = (3 * PI**2 * DT)**THD
+    KS = SQRT( 4 * KF / PI )
+    S = GDMT / (2 * KF * DT)
+    T = GDMT / (2 * KS * DT)
+    ZETA = ( D(1) - D(2) ) / DT
+    ZETA = MAX( -1.D0+DENMIN, ZETA )
+    ZETA = MIN(  1.D0-DENMIN, ZETA )
+    PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
+    F1 = ECUNIF / GAMMA / PHI**3
+    F2 = EXP(-F1)
+    A = BETA / GAMMA / (F2-1)
+    F3 = T**2 + A * T**4
+    F4 = BETA/GAMMA * F3 / (1 + A*F3)
+    F5 = 0.002568D0 + 0.023266D0*RS + 7.389D-6*RS**2
+    F6 = 1.0D0 + 8.723D0*RS + 0.472D0*RS**2 + 0.07389D0*RS**3
+    F7 = EXP(-100.0D0 * S**2 * PHI**4)
+    F8 =  15.75592D0*(0.001667212D0 + F5/F6 -0.004235D0 + &
+        3.0D0*0.001667212D0/7.0D0)
+    H = GAMMA * PHI**3 * LOG( 1 + F4 ) + F8 * T**2 * F7
+    FC = ECUNIF + H
+
+    ! Find correlation energy derivatives
+    DRSDD = - THD * RS / DT
+    DKFDD =   THD * KF / DT
+    DKSDD = HALF * KS * DKFDD / KF
+    DZDD(1) =   1 / DT - ZETA / DT
+    DZDD(2) = - 1 / DT - ZETA / DT
+    DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
+    DO IS = 1,2
+      DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
+      DPDD = DPDZ * DZDD(IS)
+      DSDD = -S*DKFDD/KF - S/DT    ! JMS: corrected, May.2014
+      DTDD = -T*DKSDD/KS - T/DT
+      DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
+      DF2DD = - F2 * DF1DD
+      DADD = - A * DF2DD / (F2-1)
+      DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
+      DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
+      DF5DD = (0.023266D0 + 2.0D0*7.389D-6*RS)*DRSDD
+      DF6DD = (8.723D0 + 2.0D0*0.472D0*RS &
+          + 3.0D0*0.07389D0*RS**2)*DRSDD
+      DF7DD = -200.0D0 * S * PHI**4 * DSDD * F7 &
+          -100.0D0 * S**2 * 4.0D0* PHI**3 * DPDD * F7
+      DF8DD = 15.75592D0 * DF5DD/F6 - 15.75592D0*F5*DF6DD / F6**2
+      DHDD = 3 * H * DPDD / PHI
+      DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
+      DHDD = DHDD + DF8DD * T**2 * F7
+      DHDD = DHDD + F8 * 2*T*DTDD *F7
+      DHDD = DHDD + F8 * T**2 * DF7DD
+
+      DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
+      DO IX = 1,3
+        DTDGD = (T / GDMT) * GDT(IX) / GDMT
+        DSDGD = (S / GDMT) * GDT(IX) / GDMT
+        DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
+        DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 
+        DF7DGD = -200.0D0 * S * PHI**4 * DSDGD * F7
+        DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
+        DHDGD = DHDGD + F8 * 2*T*DTDGD *F7 + F8 * T**2 *DF7DGD
+        DFCDGD(IX,IS) = DT * DHDGD
+      end DO
+    end DO
+
+    ! Find exchange energy and potential
+    FX = 0
+    DO IS = 1,2
+      DS(1) = MAX( DENMIN, 2 * D(IS) )
+      GDMS = MAX( GDMIN, 2 * GDM(IS) )
+      KFS = (3 * PI**2 * DS(1))**THD
+      S = GDMS / (2 * KFS * DS(1))
+      F4 = SQRT(1.0D0 + (7.7956D0*S)**2)
+      ASINHS = LOG(7.7956D0*S + F4)
+      F1 = 1.0D0 + 0.19645D0 * S * ASINHS
+      F2 = 0.2743D0 - 0.15084D0*EXP(-100.0D0*S*S)
+      F3 = 1.0D0 / (F1 + 0.004D0 * S*S*S*S)
+      F = (F1 + F2 * S*S ) * F3
+      CALL EXCHNG( IREL, 1, DS, EXUNIF, VXUNIF )
+      FX = FX + DS(1) * EXUNIF * F
+
+      DKFDD = THD * KFS / DS(1)
+      DSDD = S * ( -DKFDD/KFS - 1/DS(1) )
+      DF1DS = 0.19645D0 * ASINHS + &
+          0.19645D0 * S * 7.7956D0 / F4
+      DF2DS = 0.15084D0*200.0D0*S*EXP(-100.0D0*S*S)
+      DF3DS = - F3*F3 * (DF1DS + 4.0D0*0.004D0 * S*S*S)
+      DFDS =  DF1DS * F3 + DF2DS * S*S * F3 + 2.0D0 * S * F2 * F3 &
+          + (F1 + F2 * S*S ) * DF3DS   
+      DFXDD(IS) = VXUNIF(1) * F + DS(1) * EXUNIF * DFDS * DSDD
+
+      DO IX = 1,3
+        GDS = 2 * GD(IX,IS)
+        DSDGD = (S / GDMS) * GDS / GDMS
+        DFDGD = DFDS * DSDGD
+        DFXDGD(IX,IS) = DS(1) * EXUNIF * DFDGD
+      end DO
+    end DO
+    FX = HALF * FX / DT
+
+    ! Set output arguments
+    EX = FX
+    EC = FC
+    DO IS = 1,nspin
+      DEXDD(IS) = DFXDD(IS)
+      DECDD(IS) = DFCDD(IS)
+      DO IX = 1,3
+        DEXDGD(IX,IS) = DFXDGD(IX,IS)
+        DECDGD(IX,IS) = DFCDGD(IX,IS)
+      end DO
+    end DO
+
+  END SUBROUTINE PW91XC
+
+
+
+  subroutine blypxc(nspin,dens,gdens,EX,EC, &
+      dEXdd,dECdd,dEXdgd,dECdgd) 
+    ! ***************************************************************
+    ! Implements Becke gradient exchange functional (A.D. 
+    ! Becke, Phys. Rev. A 38, 3098 (1988)) and Lee, Yang, Parr
+    ! correlation functional (C. Lee, W. Yang, R.G. Parr, Phys. Rev. B
+    ! 37, 785 (1988)), as modificated by Miehlich,Savin,Stoll and Preuss,
+    ! Chem. Phys. Lett. 157,200 (1989). See also Johnson, Gill and Pople,
+    ! J. Chem. Phys. 98, 5612 (1993). Some errors were detected in this
+    ! last paper, so not all of the expressions correspond exactly to those
+    ! implemented here.
+    ! Written by Maider Machado. July 1998.
+    ! **************** INPUT ******************************************** 
+    ! integer nspin          : Number of spin polarizations (1 or 2)
+    ! real*8  dens(nspin)    : Total electron density (if nspin=1) or
+    !                           spin electron density (if nspin=2)
+    ! real*8  gdens(3,nspin) : Total or spin density gradient
+    ! ******** OUTPUT *****************************************************
+    ! real*8  ex             : Exchange energy density
+    ! real*8  ec             : Correlation energy density
+    ! real*8  dexdd(nspin)   : Partial derivative
+    !                           d(DensTot*Ex)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          exchange potential
+    ! real*8  decdd(nspin)   : Partial derivative
+    !                           d(DensTot*Ec)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          correlation potential
+    ! real*8  dexdgd(3,nspin): Partial derivative
+    !                           d(DensTot*Ex)/d(GradDens(i,ispin))
+    ! real*8  decdgd(3,nspin): Partial derivative
+    !                           d(DensTot*Ec)/d(GradDens(i,ispin))
+    ! ********* UNITS ****************************************************
+    ! Lengths in Bohr
+    ! Densities in electrons per Bohr**3
+    ! Energies in Hartrees
+    ! Gradient vectors in cartesian coordinates
+    ! ********************************************************************
+
+    integer nspin
+    real(dp)   dens(nspin), gdens(3,nspin), EX, EC, &
+        dEXdd(nspin), dECdd(nspin), dEXdgd(3,nspin), &
+        dECdgd(3,nspin)
+
+    ! Internal variables
+    integer is,ix
+    real(dp)   pi, beta, thd, tthd, thrhlf, half, fothd, &
+        d(2),gd(3,2),dmin, ash,gdm(2),denmin,dt, &
+        g(2),x(2),a,b,c,dd,onzthd,gdmin, &
+        ga, gb, gc,becke,dbecgd(3,2), &
+        dgdx(2), dgdxa, dgdxb, dgdxc,dgdxd,dbecdd(2), &
+        den,omega, domega, delta, ddelta,cf, &
+        gam11, gam12, gam22, LYPa, LYPb1, &
+        LYPb2,dLYP11,dLYP12,dLYP22,LYP, &
+        dd1g11,dd1g12,dd1g22,dd2g12,dd2g11,dd2g22, &
+        dLYPdd(2),dg11dd(3,2),dg22dd(3,2), &
+        dg12dd(3,2),dLYPgd(3,2)
+
+    ! Lower bounds of density and its gradient to avoid divisions by zero
+    parameter ( denmin=1.0e-8_dp )
+    parameter ( gdmin=1.0e-8_dp )
+    parameter ( dmin=1.0e-5_dp ) 
+
+    ! Fix some numerical parameters 
+    parameter ( thd = 1.d0/3.d0, tthd=2.d0/3.d0 )
+    parameter ( thrhlf=1.5d0, half=0.5d0, &
+        fothd=4.d0/3.d0, onzthd=11.d0/3.d0)
+
+    ! Empirical parameter for Becke exchange functional (a.u.)
+    parameter(beta= 0.0042d0) 
+
+    ! Constants for LYP functional (a.u.) 
+    parameter(a=0.04918d0, b=0.132d0, c=0.2533d0, dd=0.349d0)
+
+    pi= 4*atan(1.d0)
+
+    ! Translate density and its gradient to new variables
+    if (nspin .eq. 1) then
+      d(1) = half * dens(1)
+      d(1) = max(denmin,d(1))
+      d(2) = d(1)
+      dt = max( denmin, dens(1) )
+      do ix = 1,3
+        gd(ix,1) = half * gdens(ix,1)    
+        gd(ix,2) = gd(ix,1)
+      enddo
+    else
+      d(1) = dens(1)
+      d(2) = dens(2)
       do is=1,2
+        d(is) = max (denmin,d(is))
+      enddo
+      dt = max( denmin, dens(1)+dens(2) )  
+      do ix = 1,3
+        gd(ix,1) = gdens(ix,1)
+        gd(ix,2) = gdens(ix,2)
+      enddo
+    endif
+
+    gdm(1) = sqrt( gd(1,1)**2 + gd(2,1)**2 + gd(3,1)**2 )
+    gdm(2) = sqrt( gd(1,2)**2 + gd(2,2)**2 + gd(3,2)**2 )
+
+    do is=1,2
       gdm(is)= max(gdm(is),gdmin)
-      enddo
+    enddo
 
-c Find Becke exchange energy
-       ga = -thrhlf*(3.d0/4.d0/pi)**thd
-      do is=1,2
-       if(d(is).lt.dmin) then
+    ! Find Becke exchange energy
+    ga = -thrhlf*(3.d0/4.d0/pi)**thd
+    do is=1,2
+      if(d(is).lt.dmin) then
         g(is)=ga
-       else
+      else
         x(is) = gdm(is)/d(is)**fothd
         gb = beta*x(is)**2
         ash=log(x(is)+sqrt(x(is)**2+1)) 
         gc = 1+6*beta*x(is)*ash        
         g(is) = ga-gb/gc
-       endif
-      enddo
-
-c   Density of energy 
-      becke=(g(1)*d(1)**fothd+g(2)*d(2)**fothd)/dt
-
-c Exchange energy derivatives
-       do is=1,2
-        if(d(is).lt.dmin)then
-         dbecdd(is)=0.
-         do ix=1,3
+      endif
+    enddo
+
+    !   Density of energy 
+    becke=(g(1)*d(1)**fothd+g(2)*d(2)**fothd)/dt
+
+    ! Exchange energy derivatives
+    do is=1,2
+      if(d(is).lt.dmin)then
+        dbecdd(is)=0.
+        do ix=1,3
           dbecgd(ix,is)=0.
-         enddo
-        else
+        enddo
+      else
         dgdxa=6*beta**2*x(is)**2
         ash=log(x(is)+sqrt(x(is)**2+1))
         dgdxb=x(is)/sqrt(x(is)**2+1)-ash
@@ -1328,1513 +1298,1484 @@
         dgdx(is)=(dgdxa*dgdxb+dgdxc)/dgdxd
         dbecdd(is)=fothd*d(is)**thd*(g(is)-x(is)*dgdx(is))
         do ix=1,3
-         dbecgd(ix,is)=d(is)**(-fothd)*dgdx(is)*gd(ix,is)/x(is)
-        enddo 
-        endif
-       enddo
-
-c  Lee-Yang-Parr correlation energy
-      den=1+dd*dt**(-thd)
-      omega=dt**(-onzthd)*exp(-c*dt**(-thd))/den
-      delta=c*dt**(-thd)+dd*dt**(-thd)/den
-      cf=3.*(3*pi**2)**tthd/10.
-      gam11=gdm(1)**2
-      gam12=gd(1,1)*gd(1,2)+gd(2,1)*gd(2,2)+gd(3,1)*gd(3,2)
-      gam22=gdm(2)**2
-      LYPa=-4*a*d(1)*d(2)/(den*dt)
-      LYPb1=2**onzthd*cf*a*b*omega*d(1)*d(2)
-      LYPb2=d(1)**(8./3.)+d(2)**(8./3.)
-      dLYP11=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.)
-     .*d(1)/dt)-d(2)**2)
-      dLYP12=-a*b*omega*(d(1)*d(2)/9.*(47.-7.*delta)
-     .-fothd*dt**2)
-      dLYP22=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.)*
-     .d(2)/dt)-d(1)**2)
-
-c    Density of energy
-      LYP=(LYPa-LYPb1*LYPb2+dLYP11*gam11+dLYP12*gam12
-     .+dLYP22*gam22)/dt
-
-c   Correlation energy derivatives
-       domega=-thd*dt**(-fothd)*omega*(11.*dt**thd-c-dd/den)
-       ddelta=thd*(dd**2*dt**(-5./3.)/den**2-delta/dt)
-
-c   Second derivatives with respect to the density
-       dd1g11=domega/omega*dLYP11-a*b*omega*(d(2)/9.*
-     . (1.-3.*delta-2*(delta-11.)*d(1)/dt)-d(1)*d(2)/9.*
-     . ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2))
-
-       dd1g12=domega/omega*dLYP12-a*b*omega*(d(2)/9.*
-     . (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt)
-
-      dd1g22=domega/omega*dLYP22-a*b*omega*(1./9.*d(2)
-     . *(1.-3.*delta-(delta-11.)*d(2)/dt)-d(1)*d(2)/9.*
-     . ((3.+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2)-2*d(1))
-
-      dd2g22=domega/omega*dLYP22-a*b*omega*(d(1)/9.*
-     . (1.-3.*delta-2*(delta-11.)*d(2)/dt)-d(1)*d(2)/9.*
-     . ((3+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2))
-      
-      dd2g12=domega/omega*dLYP12-a*b*omega*(d(1)/9.*
-     . (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt)
-      
-      dd2g11=domega/omega*dLYP11-a*b*omega*(1./9.*d(1)
-     . *(1.-3.*delta-(delta-11.)*d(1)/dt)-d(1)*d(2)/9.*
-     . ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2)-2*d(2))
-
-        dLYPdd(1)=-4*a/den*d(1)*d(2)/dt*
-     . (thd*dd*dt**(-fothd)/den
-     . +1./d(1)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)*
-     . (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(2)*(onzthd*
-     . d(1)**(8./3.)+d(2)**(8./3.)))+dd1g11*gam11+
-     . dd1g12*gam12+dd1g22*gam22
-
-       dLYPdd(2)=-4*a/den*d(1)*d(2)/dt*(thd*dd*dt**(-fothd)/den
-     . +1./d(2)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)*
-     . (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(1)*(onzthd*
-     . d(2)**(8./3.)+d(1)**(8./3.)))+dd2g22*gam22+
-     . dd2g12*gam12+dd2g11*gam11
-
-c second derivatives with respect to the density gradient
-        do is=1,2
-          do ix=1,3
-           dg11dd(ix,is)=2*gd(ix,is)
-           dg22dd(ix,is)=2*gd(ix,is)
-          enddo
-        enddo
-        do ix=1,3
-          dLYPgd(ix,1)=dLYP11*dg11dd(ix,1)+dLYP12*gd(ix,2)
-          dLYPgd(ix,2)=dLYP22*dg22dd(ix,2)+dLYP12*gd(ix,1)
-        enddo
-
-       EX=becke
-       EC=LYP
-       do is=1,nspin
-        dEXdd(is)=dbecdd(is)
-        dECdd(is)=dLYPdd(is)
-        do ix=1,3
-         dEXdgd(ix,is)=dbecgd(ix,is)
-         dECdgd(ix,is)=dLYPgd(ix,is)
-        enddo
-       enddo
-       end subroutine blypxc
-
-
-
-      SUBROUTINE RPBEXC( IREL, nspin, Dens, GDens,
-     .                   EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Hammer's RPBE Generalized-Gradient-Approximation (GGA).
-C A revision of PBE (Perdew-Burke-Ernzerhof) 
-C Ref: Hammer, Hansen & Norskov, PRB 59, 7413 (1999) and
-C J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
-C
-C Written by M.V. Fernandez-Serra. March 2004. On the PBE routine of
-C L.C.Balbas and J.M.Soler. December 1996.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C EXCHNG, PW92C
-C ********************************************************************
-
-      implicit          none
-      INTEGER           IREL, nspin
-      real(dp)          Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
-     .                  DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
-
-C Internal variables
-      INTEGER
-     .  IS, IX
-
-      real(dp)
-     .  A, BETA, D(2), DADD, DECUDD, DENMIN, 
-     .  DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
-     .  DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
-     .  DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, 
-     .  DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), 
-     .  EC, ECUNIF, EX, EXUNIF,
-     .  F, F1, F2, F3, F4, FC, FX, FOUTHD,
-     .  GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
-     .  H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
-     .  T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      PARAMETER ( DENMIN = 1.D-12 )
-      PARAMETER ( GDMIN  = 1.D-12 )
-
-C Fix some numerical parameters
-      PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
-     .            THD=1.D0/3.D0, THRHLF=1.5D0,
-     .            TWO=2.D0, TWOTHD=2.D0/3.D0 )
-
-C Fix some more numerical constants
-      PI = 4 * ATAN(1.D0)
-      BETA = 0.066725D0
-      GAMMA = (1 - LOG(TWO)) / PI**2
-      MU = BETA * PI**2 / 3
-      KAPPA = 0.804D0
-
-C Translate density and its gradient to new variables
-      IF (nspin .EQ. 1) THEN
-        D(1) = HALF * Dens(1)
-        D(2) = D(1)
-        DT = MAX( DENMIN, Dens(1) )
-        DO 10 IX = 1,3
-          GD(IX,1) = HALF * GDens(IX,1)
-          GD(IX,2) = GD(IX,1)
-          GDT(IX) = GDens(IX,1)
-   10   CONTINUE
-      ELSE
-        D(1) = Dens(1)
-        D(2) = Dens(2)
-        DT = MAX( DENMIN, Dens(1)+Dens(2) )
-        DO 20 IX = 1,3
-          GD(IX,1) = GDens(IX,1)
-          GD(IX,2) = GDens(IX,2)
-          GDT(IX) = GDens(IX,1) + GDens(IX,2)
-   20   CONTINUE
-      ENDIF
-      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = MAX( GDMIN, GDMT )
-
-C Find local correlation energy and potential
-      CALL PW92C( 2, D, ECUNIF, VCUNIF )
-
-C Find total correlation energy
-      RS = ( 3 / (4*PI*DT) )**THD
-      KF = (3 * PI**2 * DT)**THD
-      KS = SQRT( 4 * KF / PI )
-      ZETA = ( D(1) - D(2) ) / DT
-      ZETA = MAX( -1.D0+DENMIN, ZETA )
-      ZETA = MIN(  1.D0-DENMIN, ZETA )
-      PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
-      T = GDMT / (2 * PHI * KS * DT)
-      F1 = ECUNIF / GAMMA / PHI**3
-      F2 = EXP(-F1)
-      A = BETA / GAMMA / (F2-1)
-      F3 = T**2 + A * T**4
-      F4 = BETA/GAMMA * F3 / (1 + A*F3)
-      H = GAMMA * PHI**3 * LOG( 1 + F4 )
-      FC = ECUNIF + H
-
-C Find correlation energy derivatives
-      DRSDD = - (THD * RS / DT)
-      DKFDD =   THD * KF / DT
-      DKSDD = HALF * KS * DKFDD / KF
-      DZDD(1) =   1 / DT - ZETA / DT
-      DZDD(2) = - (1 / DT) - ZETA / DT
-      DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
-      DO 40 IS = 1,2
-        DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
-        DPDD = DPDZ * DZDD(IS)
-        DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
-        DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
-        DF2DD = (- F2) * DF1DD
-        DADD = (- A) * DF2DD / (F2-1)
-        DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
-        DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
-        DHDD = 3 * H * DPDD / PHI
-        DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
-        DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
-
-        DO 30 IX = 1,3
-          DTDGD = (T / GDMT) * GDT(IX) / GDMT
-          DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
-          DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 
-          DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
-          DFCDGD(IX,IS) = DT * DHDGD
-   30   CONTINUE
-   40 CONTINUE
-
-C Find exchange energy and potential
-      FX = 0
-      DO 60 IS = 1,2
-        DS(IS)   = MAX( DENMIN, 2 * D(IS) )
-        GDMS = MAX( GDMIN, 2 * GDM(IS) )
-        KFS = (3 * PI**2 * DS(IS))**THD
-        S = GDMS / (2 * KFS * DS(IS))
-cea Hammer's RPBE (Hammer, Hansen & Norskov PRB 59 7413 (99)
-cea     F1 = DEXP( - MU * S**2 / KAPPA)
-cea     F = 1 + KAPPA * (1 - F1)
-cea Following is standard PBE
-cea     F1 = 1 + MU * S**2 / KAPPA
-cea     F = 1 + KAPPA - KAPPA / F1
-cea (If revPBE Zhang & Yang, PRL 80,890(1998),change PBE's KAPPA to 1.245)
-        F1 = DEXP( - MU * S**2 / KAPPA)
-        F = 1 + KAPPA * (1 - F1)
- 
-c       Note nspin=1 in call to exchng...
- 
-        CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
-        FX = FX + DS(IS) * EXUNIF * F
-
-cMVFS   The derivatives of F  also need to be changed for Hammer's RPBE.
-cMVFS   DF1DD = 2 * F1 * DSDD  * ( - MU * S / KAPPA)
-cMVFS   DF1DGD= 2 * F1 * DSDGD * ( - MU * S / KAPPA)
-cMVFS   DFDD  = -1 * KAPPA * DF1DD
-cMVFS   DFDGD = -1 * KAPPA * DFDGD
-
-        DKFDD = THD * KFS / DS(IS)
-        DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
-c       DF1DD = 2 * (F1-1) * DSDD / S
-c       DFDD = KAPPA * DF1DD / F1**2
-        DF1DD = 2* F1 * DSDD * ( - MU * S / KAPPA)
-        DFDD = -1 * KAPPA * DF1DD
-        DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
-
-        DO 50 IX = 1,3
-          GDS = 2 * GD(IX,IS)
-          DSDGD = (S / GDMS) * GDS / GDMS
-c         DF1DGD = 2 * MU * S * DSDGD / KAPPA
-c         DFDGD = KAPPA * DF1DGD / F1**2
-          DF1DGD =2*F1 * DSDGD * ( - MU * S / KAPPA)
-          DFDGD = -1 * KAPPA * DF1DGD
-          DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
-   50   CONTINUE
-   60 CONTINUE
-      FX = HALF * FX / DT
-
-C Set output arguments
-      EX = FX
-      EC = FC
-      DO 90 IS = 1,nspin
-        DEXDD(IS) = DFXDD(IS)
-        DECDD(IS) = DFCDD(IS)
-        DO 80 IX = 1,3
-          DEXDGD(IX,IS) = DFXDGD(IX,IS)
-          DECDGD(IX,IS) = DFCDGD(IX,IS)
-   80   CONTINUE
-   90 CONTINUE
-
-      END SUBROUTINE RPBEXC
-
-
-
-      SUBROUTINE WCXC( IREL, nspin, Dens, GDens,
-     .                  EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Wu-Cohen Generalized-Gradient-Approximation.
-C Ref: Z. Wu and R. E. Cohen PRB 73, 235116 (2006)
-C Written by Marivi Fernandez-Serra, with contributions by
-C Julian Gale and Alberto Garcia,
-C over the PBEXC subroutine of L.C.Balbas and J.M.Soler.
-C September, 2006.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C EXCHNG, PW92C
-C ********************************************************************
-
-      implicit          none
-      INTEGER           IREL, nspin
-      real(dp)          Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
-     .                  DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
-
-C Internal variables
-      INTEGER
-     .  IS, IX
-
-      real(dp)
-     .  A, BETA, D(2), DADD, DECUDD, DENMIN, 
-     .  DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
-     .  DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
-     .  DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, 
-     .  DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), 
-     .  XWC, DXWCDS, CWC,
-     .  EC, ECUNIF, EX, EXUNIF,
-     .  F, F1, F2, F3, F4, FC, FX, FOUTHD,
-     .  GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
-     .  H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
-     .  TEN81, 
-     .  T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      PARAMETER ( DENMIN = 1.D-12 )
-      PARAMETER ( GDMIN  = 1.D-12 )
-
-C Fix some numerical parameters
-      PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
-     .            THD=1.D0/3.D0, THRHLF=1.5D0,
-     .            TWO=2.D0, TWOTHD=2.D0/3.D0 )
-      PARAMETER ( TEN81 = 10.0d0/81.0d0 )
-
-C Fix some more numerical constants
-      PI = 4 * ATAN(1.D0)
-      BETA = 0.066725D0
-      GAMMA = (1 - LOG(TWO)) / PI**2
-      MU = BETA * PI**2 / 3
-      KAPPA = 0.804D0
-      CWC = 0.0079325D0
-
-C Translate density and its gradient to new variables
-      IF (nspin .EQ. 1) THEN
-        D(1) = HALF * Dens(1)
-        D(2) = D(1)
-        DT = MAX( DENMIN, Dens(1) )
-        DO 10 IX = 1,3
-          GD(IX,1) = HALF * GDens(IX,1)
-          GD(IX,2) = GD(IX,1)
-          GDT(IX) = GDens(IX,1)
-   10   CONTINUE
-      ELSE
-        D(1) = Dens(1)
-        D(2) = Dens(2)
-        DT = MAX( DENMIN, Dens(1)+Dens(2) )
-        DO 20 IX = 1,3
-          GD(IX,1) = GDens(IX,1)
-          GD(IX,2) = GDens(IX,2)
-          GDT(IX) = GDens(IX,1) + GDens(IX,2)
-   20   CONTINUE
-      ENDIF
-      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = MAX( GDMIN, GDMT )
-
-C Find local correlation energy and potential
-      CALL PW92C( 2, D, ECUNIF, VCUNIF )
-
-C Find total correlation energy
-      RS = ( 3 / (4*PI*DT) )**THD
-      KF = (3 * PI**2 * DT)**THD
-      KS = SQRT( 4 * KF / PI )
-      ZETA = ( D(1) - D(2) ) / DT
-      ZETA = MAX( -1.D0+DENMIN, ZETA )
-      ZETA = MIN(  1.D0-DENMIN, ZETA )
-      PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
-      T = GDMT / (2 * PHI * KS * DT)
-      F1 = ECUNIF / GAMMA / PHI**3
-      F2 = EXP(-F1)
-      A = BETA / GAMMA / (F2-1)
-      F3 = T**2 + A * T**4
-      F4 = BETA/GAMMA * F3 / (1 + A*F3)
-      H = GAMMA * PHI**3 * LOG( 1 + F4 )
-      FC = ECUNIF + H
-
-C Find correlation energy derivatives
-      DRSDD = - (THD * RS / DT)
-      DKFDD =   THD * KF / DT
-      DKSDD = HALF * KS * DKFDD / KF
-      DZDD(1) =   1 / DT - ZETA / DT
-      DZDD(2) = - (1 / DT) - ZETA / DT
-      DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
-      DO 40 IS = 1,2
-        DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
-        DPDD = DPDZ * DZDD(IS)
-        DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
-        DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
-        DF2DD = (- F2) * DF1DD
-        DADD = (- A) * DF2DD / (F2-1)
-        DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
-        DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
-        DHDD = 3 * H * DPDD / PHI
-        DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
-        DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
-
-        DO 30 IX = 1,3
-          DTDGD = (T / GDMT) * GDT(IX) / GDMT
-          DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
-          DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 
-          DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
-          DFCDGD(IX,IS) = DT * DHDGD
-   30   CONTINUE
-   40 CONTINUE
-
-C Find exchange energy and potential
-      FX = 0
-      DO 60 IS = 1,2
-        DS(IS)   = MAX( DENMIN, 2 * D(IS) )
-        GDMS = MAX( GDMIN, 2 * GDM(IS) )
-        KFS = (3 * PI**2 * DS(IS))**THD
-        S = GDMS / (2 * KFS * DS(IS))
-c
-c For PBE: 
-c
-c       x = MU * S**2
-c       dxds = 2*MU*S
-c
-c Wu-Cohen form:
-c
-        XWC= TEN81 * s**2 + (MU- TEN81) * 
-     .       S**2 * exp(-S**2) + log(1+ CWC * S**4)
-        DXWCDS = 2 * TEN81 * S + (MU - TEN81) * exp(-S**2) *
-     .           2*S * (1 - S*S) + 4 * CWC * S**3 / (1 + CWC * S**4) 
-c-------------------
-
-        F1 = 1 +  XWC / KAPPA
-        F = 1 + KAPPA - KAPPA / F1
-c
-c       Note nspin=1 in call to exchng...
-c
-        CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
-        FX = FX + DS(IS) * EXUNIF * F
-
-        DKFDD = THD * KFS / DS(IS)
-        DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
-        DF1DD = DXWCDS * DSDD / KAPPA 
-        DFDD = KAPPA * DF1DD / F1**2
-        DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
-
-        DO 50 IX = 1,3
-          GDS = 2 * GD(IX,IS)
-          DSDGD = (S / GDMS) * GDS / GDMS
-          DF1DGD = DXWCDS * DSDGD / KAPPA
-          DFDGD = KAPPA * DF1DGD / F1**2
-          DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
-   50   CONTINUE
-   60 CONTINUE
-      FX = HALF * FX / DT
-
-C Set output arguments
-      EX = FX
-      EC = FC
-      DO 90 IS = 1,nspin
-        DEXDD(IS) = DFXDD(IS)
-        DECDD(IS) = DFCDD(IS)
-        DO 80 IX = 1,3
-          DEXDGD(IX,IS) = DFXDGD(IX,IS)
-          DECDGD(IX,IS) = DFCDGD(IX,IS)
-   80   CONTINUE
-   90 CONTINUE
-
-      END SUBROUTINE WCXC
-
-
-
-      SUBROUTINE AM05XC( IREL, nspin, Dens, GDens,
-     .                   EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements the Armiento Mattsson AM05 GGA.
-C Ref: R. Armiento and A. E. Mattsson, PRB 72, 085108 (2005)
-C Written by L.C.Balbas and J.M.Soler originally for PBE. December 1996. 
-C Modified by J.D. Gale for AM05. May 2009.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C am05wbs
-C ********************************************************************
-
-      use precision, only : dp
-      use am05,      only : am05wbs
-
-      implicit          none
-      integer           irel, nspin
-      real(dp)          Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
-     .                  DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
-
-C Internal variables
-      integer
-     .  is, ix
-
-      real(dp)
-     .  D(2), DENMIN, DFXDD(2), DFCDD(2), DFCDGD(3,2), 
-     .  DFXDGD(3,2), DFXDG(2), DFCDG(2),
-     .  DS(2), DT, EC, EX, FX, FC,
-     .  GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3)
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      parameter ( DENMIN = 1.D-12 )
-      parameter ( GDMIN  = 1.D-12 )
-
-C Translate density and its gradient to new variables
-      if (nspin .eq. 1) then
-        D(1) = 0.5_dp*Dens(1)
-        D(2) = D(1)
-        DT = max( DENMIN, Dens(1) )
-        do ix = 1,3
-          GD(ix,1) = 0.5_dp*GDens(ix,1)
-          GD(ix,2) = GD(ix,1)
-          GDT(ix) = GDens(ix,1)
-        enddo
-      else
-        D(1) = Dens(1)
-        D(2) = Dens(2)
-        DT = max( DENMIN, Dens(1)+Dens(2) )
-        do ix = 1,3
-          GD(ix,1) = GDens(ix,1)
-          GD(ix,2) = GDens(ix,2)
-          GDT(ix) = GDens(ix,1) + GDens(ix,2)
+          dbecgd(ix,is)=d(is)**(-fothd)*dgdx(is)*gd(ix,is)/x(is)
         enddo
       endif
-      GDM(1) = sqrt( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = sqrt( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = sqrt( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = max( GDMIN, GDMT )
-
-      D(1) = max(D(1),denmin)
-      D(2) = max(D(2),denmin)
-
-C Call AM05 subroutine
-      call am05wbs(D(1), D(2), GDM(1), GDM(2), FX, FC,
-     .             DFXDD(1), DFXDD(2), DFCDD(1), DFCDD(2), 
-     .             DFXDG(1), DFXDG(2), DFCDG(1), DFCDG(2))
-
-C Convert gradient terms into vectors
-      do is = 1,nspin
-        do ix = 1,3
-          DFXDGD(ix,is) = DFXDG(is)*GD(ix,is)
-          DFCDGD(ix,is) = DFCDG(is)*GD(ix,is)
-        enddo
-      enddo
-
-C Convert FX to form required by SIESTA - note factor of 1/2
-C is already applied in am05 code.
-      FX = FX / DT
-
-C Set output arguments
-      EX = FX
-      EC = FC
-      do is = 1,nspin
-        DEXDD(is) = DFXDD(is)
-        DECDD(is) = DFCDD(is)
-        do ix = 1,3
-          DEXDGD(ix,is) = DFXDGD(ix,is)
-          DECDGD(ix,is) = DFCDGD(ix,is)
-        enddo
-      enddo
-
-      END SUBROUTINE AM05XC
-
-
-
-      SUBROUTINE PW86formX( a, b, c, iRel, nSpin, Dens, GDens,
-     .                      EX, dEXdD, dEXdGD )
-
-C *********************************************************************
-C Implements Perdew-Wang-86 Generalized-Gradient-Approximation exchange-only
-C functional form, eps_x=eps_x_LDA*(1+15*a*s**2+b*s**4+c*s**6)**(1/15),
-C but with variable values for parameters a, b, and c
-C Refs: J.P.Perdew & Y.Wang, PRB 33, 8800 (1986)
-C       E.D.Murray, K.Lee & D.C.Langreth, JCTC 5, 2754 (2009)
-C Written by J.M.Soler. March 2010.
-C ******** INPUT ******************************************************
-C REAL*8  a, b, c        : Parameter a, b, and c of the PW86 functional
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C EXCHNG
-C ********************************************************************
-
-      implicit none
-
-C Input
-      real(dp),intent(in) :: a              ! Parameter of PW86 functional
-      real(dp),intent(in) :: b              ! Parameter of PW86 functional
-      real(dp),intent(in) :: c              ! Parameter of PW86 functional
-      integer, intent(in) :: iRel           ! Relativistic exchange? 0=no, 1=yes
-      integer, intent(in) :: nSpin          ! Number of spin components
-      real(dp),intent(in) :: Dens(nSpin)    ! Local electron (spin) density
-      real(dp),intent(in) :: GDens(3,nSpin) ! Gradient of electron density
-
-C Output
-      real(dp),intent(out):: EX             ! Exchange energy per electron
-      real(dp),intent(out):: dEXdD(nSpin)   ! dEx/dDens, Ex=Int(dens*epsX)
-      real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens)
-
-C Internal variables
-      INTEGER
-     .  IS, IX
-      real(dp)
-     .  D(2), DENMIN, DF1DS, DFDD, DFDGD, DFDS, DFXDD(2), DFXDGD(3,2), 
-     .  DKFDD, DS(2), DSDD, DSDGD, DT, ECUNIF, EXUNIF, F, F1, FX,
-     .  GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
-     .  KFS, PI, S, VXUNIF(2), ZETA
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      PARAMETER ( DENMIN = 1.D-12 )
-      PARAMETER ( GDMIN  = 1.D-12 )
-
-C Fix some more numerical constants
-      PI = 4 * ATAN(1.D0)
-
-C Translate density and its gradient to new variables
-      IF (nspin .EQ. 1) THEN
-        D(1) = Dens(1) / 2
-        D(2) = D(1)
-        DT = MAX( DENMIN, Dens(1) )
-        DO 10 IX = 1,3
-          GD(IX,1) = GDens(IX,1) / 2
-          GD(IX,2) = GD(IX,1)
-          GDT(IX) = GDens(IX,1)
-   10   CONTINUE
-      ELSE
-        D(1) = Dens(1)
-        D(2) = Dens(2)
-        DT = MAX( DENMIN, Dens(1)+Dens(2) )
-        DO 20 IX = 1,3
-          GD(IX,1) = GDens(IX,1)
-          GD(IX,2) = GDens(IX,2)
-          GDT(IX) = GDens(IX,1) + GDens(IX,2)
-   20   CONTINUE
-      ENDIF
-      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = MAX( GDMIN, GDMT )
-
-C Find exchange energy and potential
-      FX = 0
-      DO 60 IS = 1,2
-        DS(IS)   = MAX( DENMIN, 2 * D(IS) )
-        GDMS = MAX( GDMIN, 2 * GDM(IS) )
-        KFS = (3 * PI**2 * DS(IS))**(1._dp/3)
-        S = GDMS / (2 * KFS * DS(IS))
-        F1 = 1 + 15*a*S**2 + b*S**4 + c*S**6
-        F = F1**(1._dp/15)
-c
-c       Note nspin=1 in call to exchng...
-c
-        CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
-        FX = FX + DS(IS) * EXUNIF * F
-
-        DKFDD = KFS / DS(IS) / 3
-        DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
-        DF1DS = 30*a*S + 4*b*S**3 + 6*c*S**5
-        DFDS = F/F1/15 * DF1DS
-        DFDD = DFDS * DSDD
-        DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
-
-        DO 50 IX = 1,3
-          GDS = 2 * GD(IX,IS)
-          DSDGD = (S / GDMS) * GDS / GDMS
-          DFDGD = DFDS * DSDGD
-          DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
-   50   CONTINUE
-   60 CONTINUE
-      FX = FX / DT / 2
-
-C Set output arguments
-      EX = FX
-      DO 90 IS = 1,nspin
-        DEXDD(IS) = DFXDD(IS)
-        DO 80 IX = 1,3
-          DEXDGD(IX,IS) = DFXDGD(IX,IS)
-   80   CONTINUE
-   90 CONTINUE
-
-      END SUBROUTINE PW86formX
-
-
-
-      SUBROUTINE PW86X( IREL, nspin, Dens, GDens,
-     .                  EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Wang-86 Generalized-Gradient-Approximation 
-C exchange-only functional. Correlation energy returns as zero.
-C Ref: J.P.Perdew & Y.Wang, PRB 33, 8800 (1986)
-C Written by J.M.Soler. March 2010.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C PW86formX
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal parameters of the PW86 exchange functional
-      real(dp),parameter:: a = 0.0864_dp
-      real(dp),parameter:: b = 14.0_dp
-      real(dp),parameter:: c = 0.2_dp
-
-C Call PW86 routine with appropriate values for a, b, and c.
-      CALL PW86formX( a, b, c, IREL, nspin, Dens, GDens,
-     .                EX, DEXDD, DEXDGD )
-
-C Set correlation energy and derivatives to zero
-      EC = 0
-      DECDD(:) = 0
-      DECDGD(:,:) = 0
-
-      END SUBROUTINE PW86X
-
-
-
-      SUBROUTINE PW86RX( IREL, nspin, Dens, GDens,
-     .                  EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Perdew-Wang-86 Generalized-Gradient-Approximation 
-C exchange-only functional with the refitted parameters of
-C Murray, Lee, and Langreth. Correlation energy returns as zero.
-C Refs: J.P.Perdew & Y.Wang, PRB 33, 8800 (1986)
-C       E.D.Murray, K.Lee & D.C.Langreth, JCTC 5, 2754 (2009)
-C Written by J.M.Soler. March 2010.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C PW86formX
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal parameters of the refitted PW86 exchange functional
-      real(dp),parameter:: a = 0.1234_dp
-      real(dp),parameter:: b = 17.33_dp
-      real(dp),parameter:: c = 0.163_dp
-
-C Call PW86 routine with appropriate values for a, b, and c.
-      CALL PW86formX( a, b, c, IREL, nspin, Dens, GDens,
-     .                EX, DEXDD, DEXDGD )
-
-C Set correlation energy and derivatives to zero
-      EC = 0
-      DECDD(:) = 0
-      DECDGD(:,:) = 0
-
-      END SUBROUTINE PW86RX
-
-
-      SUBROUTINE BHX( IREL, nspin, Dens, GDens,
-     .                EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements the combination by Berland and Hyldgaard of Perdew-Wang-86 
-C Generalized-Gradient-Approximation exchange-only functional with the 
-C refitted parameters of Murray, Lee, and Langreth and Langreth-Vosko 
-C screened exchange. Correlation energy returns as zero.
-C Refs: J.P.Perdew & Y.Wang, PRB 33, 8800 (1986)
-C       E.D.Murray, K.Lee & D.C.Langreth, JCTC 5, 2754 (2009)
-C       K.Berland & P.Hyldgaard, PRB 89, 035412 (2014)
-C Written by Michelle Fritz Feb. 2014.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER NSPIN          : Number of spin polarizations (1 or 2)
-C REAL*8  DENS(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C EXCHNG
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables
-      INTEGER
-     .  IS, IX
-      real(dp)
-     .  D(2), DENMIN, DF1DS, DFDD, DFDGD, DFDS, DFLVDS, DFPW86RDS, 
-     .  DFXDD(2), DFXDGD(3,2), DKFDD, DS(2), DSDD, DSDGD, DT, ECUNIF, 
-     .  EXUNIF, FPW86R, FLV, F, F1, FX, GD(3,2), GDM(2), GDMIN, GDMS, 
-     .  GDMT, GDS, GDT(3), KFS, PI, S, VXUNIF(2), ZETA, MLV
-
-C Internal parameters of the refitted PW86 exchange functional
-      real(dp),parameter:: a = 0.1234_dp
-      real(dp),parameter:: b = 17.33_dp
-      real(dp),parameter:: c = 0.163_dp
-      real(dp),parameter:: zab = -0.8491_dp
-      real(dp),parameter:: alpha = 0.02178_dp
-      real(dp),parameter:: beta = 1.15_dp
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      PARAMETER ( DENMIN = 1.D-12 )
-      PARAMETER ( GDMIN  = 1.D-12 )
-
-C Fix some more numerical constants
-      PI = 4 * ATAN(1.D0)
-
-C Translate density and its gradient to new variables
-      IF (NSPIN .EQ. 1) THEN
-        D(1) = DENS(1) / 2
-        D(2) = D(1)
-        DT = MAX( DENMIN, DENS(1) )
-        DO 10 IX = 1,3
-          GD(IX,1) = GDENS(IX,1) / 2
-          GD(IX,2) = GD(IX,1)
-          GDT(IX) = GDENS(IX,1)
-   10   CONTINUE
-      ELSE
-        D(1) = DENS(1)
-        D(2) = DENS(2)
-        DT = MAX( DENMIN, DENS(1)+DENS(2) )
-        DO 20 IX = 1,3
-          GD(IX,1) = GDENS(IX,1)
-          GD(IX,2) = GDENS(IX,2)
-          GDT(IX) = GDENS(IX,1) + GDENS(IX,2)
-   20   CONTINUE
-      ENDIF
-      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = MAX( GDMIN, GDMT )
-
-C Find exchange energy and potential
-      FX = 0
-      DO 60 IS = 1,2
-        DS(IS)   = MAX( DENMIN, 2 * D(IS) )
-        GDMS = MAX( GDMIN, 2 * GDM(IS) )
-        KFS = (3 * PI**2 * DS(IS))**(1._dp/3)
-        S = GDMS / (2 * KFS * DS(IS))
-        F1 = 1 + 15*a*S**2 + b*S**4 + c*S**6
-        FPW86R = F1**(1._dp/15)
-        MLV = -zab/9._dp
-        FLV = 1 + MLV*S**2
-        F = (1 / (1 + alpha*S**6)) * FLV
-        F = F + ((alpha*S**6) / (beta + alpha*S**6)) * FPW86R
-        
-c
-c       Note nspin=1 in call to exchng...
-c
-        CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
-        FX = FX + DS(IS) * EXUNIF * F
-
-        DKFDD = KFS / DS(IS) / 3
-        DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
-        DF1DS = 30*a*S + 4*b*S**3 + 6*c*S**5
-        DFLVDS = 2*MLV*S
-        DFPW86RDS = FPW86R/F1/15 * DF1DS
-        DFDS = -((6*alpha*S**5) / (1 + alpha*S**6)**2) * FLV
-        DFDS = DFDS + (1 / (1 + alpha*S**6)) * DFLVDS
-        DFDS = DFDS + ((6*alpha*S**5) / (beta + alpha*S**6)) * FPW86R
-        DFDS =DFDS-((6*alpha**2*S**11)/(beta + alpha*S**6)**2)*FPW86R
-        DFDS = DFDS + ((alpha*S**6) / (beta + alpha*S**6)) * DFPW86RDS
-        DFDD = DFDS * DSDD
-        DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
-
-        DO 50 IX = 1,3
-          GDS = 2 * GD(IX,IS)
-          DSDGD = (S / GDMS) * GDS / GDMS
-          DFDGD = DFDS * DSDGD
-          DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
-   50   CONTINUE
-   60 CONTINUE
-      FX = FX / DT / 2
-
-      EX = FX
-      DO 90 IS = 1,NSPIN
-        DEXDD(IS) = DFXDD(IS)
-        DO 80 IX = 1,3
-          DEXDGD(IX,IS) = DFXDGD(IX,IS)
-   80   CONTINUE
-   90 CONTINUE
-
-C Set correlation energy and derivatives to zero
-      EC = 0
-      DECDD(:) = 0
-      DECDGD(:,:) = 0
-
-      END SUBROUTINE BHX
-
-
-      SUBROUTINE B88formX( beta, mu, c, iRel, nSpin, Dens, GDens,
-     .                     EX, dEXdD, dEXdGD )
-
-C *********************************************************************
-C Implements Becke-88 Generalized-Gradient-Approximation exchange-only
-C functional form, eps_x=eps_x_LDA*(1+mu*s**2/(1+beta*s*asinh(c*s))),
-C but with variable values for parameters beta, mu, and c
-C Refs: A.D.Becke, PRA 38, 3098 (1988)
-C       J.Klimes, D.R.Bowler, and A.Michaelides, JPCM 22, 022201 (2009)
-C Written by J.M.Soler. April 2010.
-C ******** INPUT ******************************************************
-C REAL*8  beta, mu, c    : Parameter beta mu, and c of the B88 functional,
-C                          as formulated by Klimes et al
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C EXCHNG
-C ********************************************************************
-
-      implicit none
-
-C Input
-      real(dp),intent(in) :: beta           ! Parameter of B88 functional
-      real(dp),intent(in) :: mu             ! Parameter of B88 functional
-      real(dp),intent(in) :: c              ! Parameter of B88 functional
-      integer, intent(in) :: iRel           ! Relativistic exchange? 0=no, 1=yes
-      integer, intent(in) :: nSpin          ! Number of spin components
-      real(dp),intent(in) :: Dens(nSpin)    ! Local electron (spin) density
-      real(dp),intent(in) :: GDens(3,nSpin) ! Gradient of electron density
-
-C Output
-      real(dp),intent(out):: EX             ! Exchange energy per electron
-      real(dp),intent(out):: dEXdD(nSpin)   ! dEx/dDens, Ex=Int(dens*epsX)
-      real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens)
-
-C Internal variables
-      INTEGER
-     .  IS, IX
-      real(dp)
-     .  ASINHCS, CS, D(2), DASINHDS, 
-     .  DENMIN, DF1DS, DFDD, DFDGD, DFDS, DFXDD(2), DFXDGD(3,2), 
-     .  DKFDD, DS(2), DSDD, DSDGD, DT, ECUNIF, EXUNIF, F, F1, FX,
-     .  GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
-     .  KFS, PI, S, VXUNIF(2), ZETA
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      PARAMETER ( DENMIN = 1.D-12 )
-      PARAMETER ( GDMIN  = 1.D-12 )
-
-C Fix some more numerical constants
-      PI = 4 * ATAN(1.D0)
-
-C Translate density and its gradient to new variables
-      IF (nspin .EQ. 1) THEN
-        D(1) = Dens(1) / 2
-        D(2) = D(1)
-        DT = MAX( DENMIN, Dens(1) )
-        DO 10 IX = 1,3
-          GD(IX,1) = GDens(IX,1) / 2
-          GD(IX,2) = GD(IX,1)
-          GDT(IX) = GDens(IX,1)
-   10   CONTINUE
-      ELSE
-        D(1) = Dens(1)
-        D(2) = Dens(2)
-        DT = MAX( DENMIN, Dens(1)+Dens(2) )
-        DO 20 IX = 1,3
-          GD(IX,1) = GDens(IX,1)
-          GD(IX,2) = GDens(IX,2)
-          GDT(IX) = GDens(IX,1) + GDens(IX,2)
-   20   CONTINUE
-      ENDIF
-      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = MAX( GDMIN, GDMT )
-
-C Find exchange energy and potential
-      FX = 0
-      DO 60 IS = 1,2
-        DS(IS)   = MAX( DENMIN, 2 * D(IS) )
-        GDMS = MAX( GDMIN, 2 * GDM(IS) )
-        KFS = (3 * PI**2 * DS(IS))**(1._dp/3)
-        S = GDMS / (2 * KFS * DS(IS))
-        CS = C * S
-        ! Next line is ASINH(CS). See Numerical Recipes
-        ASINHCS = log(CS+sqrt(CS**2+1))
-        F1 = 1 + beta * S * ASINHCS
-        F = 1 + mu * S**2 / F1
-c
-c       Note nspin=1 in call to exchng...
-c
-        CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
-        FX = FX + DS(IS) * EXUNIF * F
-
-        DKFDD = KFS / DS(IS) / 3
-        DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
-        DASINHDS = C * (1 + CS/sqrt(CS**2+1)) / (CS+sqrt(CS**2+1))
-        DF1DS = beta * ASINHCS + beta * S * DASINHDS
-        DFDS = 2*mu*S/F1 - mu*S**2/F1**2 * DF1DS
-        DFDD = DFDS * DSDD
-        DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
-
-        DO 50 IX = 1,3
-          GDS = 2 * GD(IX,IS)
-          DSDGD = (S / GDMS) * GDS / GDMS
-          DFDGD = DFDS * DSDGD
-          DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
-   50   CONTINUE
-   60 CONTINUE
-      FX = FX / DT / 2
-
-C Set output arguments
-      EX = FX
-      DO 90 IS = 1,nspin
-        DEXDD(IS) = DFXDD(IS)
-        DO 80 IX = 1,3
-          DEXDGD(IX,IS) = DFXDGD(IX,IS)
-   80   CONTINUE
-   90 CONTINUE
-
-      END SUBROUTINE B88formX
-
-
-
-      SUBROUTINE B88X( IREL, nspin, Dens, GDens,
-     .                 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Becke-88 Generalized-Gradient-Approximation 
-C exchange-only functional. Correlation energy returns as zero.
-C Refs: A.D.Becke, PRA 38, 3098 (1988)
-C       J.Klimes, D.R.Bowler, and A.Michaelides, JPCM 22, 022201 (2009)
-C Written by J.M.Soler. April 2010.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C PW86formX
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables and arrays
-      real(dp):: beta, c, mu, pi
-
-C Internal parameters of the B88 exchange functional, written as
-C by Klimes et al
-      pi = acos(-1._dp)
-      c = 2*(6*pi**2)**(1._dp/3)
-      mu = 0.2743_dp
-      beta = 9 * mu * (6/pi)**(1._dp/3) / (2*c)
-
-C Call PW86 routine with appropriate values for a, b, and c.
-      CALL B88formX( beta, mu, c, IREL, nspin, Dens, GDens,
-     .               EX, DEXDD, DEXDGD )
-
-C Set correlation energy and derivatives to zero
-      EC = 0
-      DECDD(:) = 0
-      DECDGD(:,:) = 0
-
-      END SUBROUTINE B88X
-
-
-
-      SUBROUTINE B88KBMX( IREL, nspin, Dens, GDens,
-     .                    EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Becke-88 GGA exchange functional, reparametrized by Klimes
-C et al for its combined use with vdW-DF nonlocal correlation.
-C Correlation energy returns as zero.
-C Refs: A.D.Becke, PRA 38, 3098 (1988)
-C       J.Klimes, D.R.Bowler, and A.Michaelides, JPCM 22, 022201 (2009)
-C Written by J.M.Soler. April 2010.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C PW86formX
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables and arrays
-      real(dp):: beta, c, mu, pi
-
-C Klimes et al parameters for the B88 exchange functional
-      pi = acos(-1._dp)
-      c = 2*(6*pi**2)**(1._dp/3)
-      mu = 0.22_dp
-      beta = mu / 1.2_dp
-
-C Call PW86 routine with appropriate values for a, b, and c.
-      CALL B88formX( beta, mu, c, IREL, nspin, Dens, GDens,
-     .               EX, DEXDD, DEXDGD )
-
-C Set correlation energy and derivatives to zero
-      EC = 0
-      DECDD(:) = 0
-      DECDGD(:,:) = 0
-
-      END SUBROUTINE B88KBMX
-
-
-      SUBROUTINE C09X( IREL, nspin, Dens, GDens,
-     .                 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
-
-C *********************************************************************
-C Implements Cooper-09 exchange-only GGA (for use with vdW-DF1 correlation). 
-C Correlation energy returns as zero.
-C Ref: V.R.Cooper, PRB 81, 161104(R) (2010)
-C Written by J.M.Soler. April 2010.
-C ******** INPUT ******************************************************
-C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
-C INTEGER nspin          : Number of spin polarizations (1 or 2)
-C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
-C                           spin electron density (if nspin=2)
-C REAL*8  GDens(3,nspin) : Total or spin density gradient
-C ******** OUTPUT *****************************************************
-C REAL*8  EX             : Exchange energy density
-C REAL*8  EC             : Correlation energy density
-C REAL*8  DEXDD(nspin)   : Partial derivative
-C                           d(DensTot*Ex)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          exchange potential
-C REAL*8  DECDD(nspin)   : Partial derivative
-C                           d(DensTot*Ec)/dDens(ispin),
-C                           where DensTot = Sum_ispin( Dens(ispin) )
-C                          For a constant density, this is the
-C                          correlation potential
-C REAL*8  DEXDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ex)/d(GradDens(i,ispin))
-C REAL*8  DECDGD(3,nspin): Partial derivative
-C                           d(DensTot*Ec)/d(GradDens(i,ispin))
-C ********* UNITS ****************************************************
-C Lengths in Bohr
-C Densities in electrons per Bohr**3
-C Energies in Hartrees
-C Gradient vectors in cartesian coordinates
-C ********* ROUTINES CALLED ******************************************
-C none
-C ********************************************************************
-
-      IMPLICIT NONE
-
-C Passed arguments
-      integer, intent(in) :: IREL, NSPIN
-      real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN)
-      real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN),
-     .                       DEXDD(NSPIN), DEXDGD(3,NSPIN)
-
-C Internal variables and arrays (V.R.Cooper, PRB 81, 161104(R) (2010))
-      real(dp),parameter:: alpha = 0.0483_dp
-      real(dp),parameter:: kappa = 1.245_dp
-      real(dp),parameter:: mu    = 0.0617
-
-C Internal variables
-      INTEGER
-     .  IS, IX
-      real(dp)
-     .  D(2), DENMIN, DES2DS, DF1DS, DFDD, DFDGD, DFDS,  
-     .  DFXDD(2), DFXDGD(3,2), DKFDD, DS(2), DSDD, DSDGD, DT, 
-     .  ECUNIF, ES2, EXUNIF, F, FX,
-     .  GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
-     .  KFS, PI, S, S2, VXUNIF(2)
-
-C Lower bounds of density and its gradient to avoid divisions by zero
-      PARAMETER ( DENMIN = 1.D-12 )
-      PARAMETER ( GDMIN  = 1.D-12 )
-
-C Fix some more numerical constants
-      PI = 4 * ATAN(1.D0)
-
-C Translate density and its gradient to new variables
-      IF (nspin .EQ. 1) THEN
-        D(1) = Dens(1) / 2
-        D(2) = D(1)
-        DT = MAX( DENMIN, Dens(1) )
-        DO 10 IX = 1,3
-          GD(IX,1) = GDens(IX,1) / 2
-          GD(IX,2) = GD(IX,1)
-          GDT(IX) = GDens(IX,1)
-   10   CONTINUE
-      ELSE
-        D(1) = Dens(1)
-        D(2) = Dens(2)
-        DT = MAX( DENMIN, Dens(1)+Dens(2) )
-        DO 20 IX = 1,3
-          GD(IX,1) = GDens(IX,1)
-          GD(IX,2) = GDens(IX,2)
-          GDT(IX) = GDens(IX,1) + GDens(IX,2)
-   20   CONTINUE
-      ENDIF
-      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
-      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
-      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
-      GDMT = MAX( GDMIN, GDMT )
-
-C Find exchange energy and potential
-      FX = 0
-      DO 60 IS = 1,2
-        DS(IS) = MAX( DENMIN, 2 * D(IS) )
-        GDMS = MAX( GDMIN, 2 * GDM(IS) )
-        KFS = (3 * PI**2 * DS(IS))**(1._dp/3)
-        S = GDMS / (2 * KFS * DS(IS))
-        S2 = S**2
-        
-        ! Next lines are the core of the C09 functional
-        ES2 = exp(-alpha*S2/2)
-        F = 1 + mu*S2*ES2**2 + kappa*(1-ES2)
-        DES2DS = -alpha*S*ES2
-        DFDS = 2*mu*S*ES2**2 + 2*mu*S2*ES2*DES2DS - kappa*DES2DS
-
-        ! Note nspin=1 in call to exchng...
-        CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
-        FX = FX + DS(IS) * EXUNIF * F
-
-        DKFDD = KFS / DS(IS) / 3
-        DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
-        DFDD = DFDS * DSDD
-        DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
-
-        DO 50 IX = 1,3
-          GDS = 2 * GD(IX,IS)
-          DSDGD = (S / GDMS) * GDS / GDMS
-          DFDGD = DFDS * DSDGD
-          DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
-   50   CONTINUE
-   60 CONTINUE
-      FX = FX / DT / 2
-
-C Set output arguments
-      EX = FX
-      DO 90 IS = 1,nspin
-        DEXDD(IS) = DFXDD(IS)
-        DO 80 IX = 1,3
-          DEXDGD(IX,IS) = DFXDGD(IX,IS)
-   80   CONTINUE
-   90 CONTINUE
-
-C Set correlation energy and derivatives to zero
-      EC = 0
-      DECDD(:) = 0
-      DECDGD(:,:) = 0
-
-      END SUBROUTINE C09X
-
-      END MODULE m_ggaxc
+    enddo
+
+    !  Lee-Yang-Parr correlation energy
+    den=1+dd*dt**(-thd)
+    omega=dt**(-onzthd)*exp(-c*dt**(-thd))/den
+    delta=c*dt**(-thd)+dd*dt**(-thd)/den
+    cf=3.*(3*pi**2)**tthd/10.
+    gam11=gdm(1)**2
+    gam12=gd(1,1)*gd(1,2)+gd(2,1)*gd(2,2)+gd(3,1)*gd(3,2)
+    gam22=gdm(2)**2
+    LYPa=-4*a*d(1)*d(2)/(den*dt)
+    LYPb1=2**onzthd*cf*a*b*omega*d(1)*d(2)
+    LYPb2=d(1)**(8./3.)+d(2)**(8./3.)
+    dLYP11=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.) &
+        *d(1)/dt)-d(2)**2)
+    dLYP12=-a*b*omega*(d(1)*d(2)/9.*(47.-7.*delta) &
+        -fothd*dt**2)
+    dLYP22=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.)* &
+        d(2)/dt)-d(1)**2)
+
+    !    Density of energy
+    LYP=(LYPa-LYPb1*LYPb2+dLYP11*gam11+dLYP12*gam12 &
+        +dLYP22*gam22)/dt
+
+    !   Correlation energy derivatives
+    domega=-thd*dt**(-fothd)*omega*(11.*dt**thd-c-dd/den)
+    ddelta=thd*(dd**2*dt**(-5./3.)/den**2-delta/dt)
+
+    !   Second derivatives with respect to the density
+    dd1g11=domega/omega*dLYP11-a*b*omega*(d(2)/9.* &
+        (1.-3.*delta-2*(delta-11.)*d(1)/dt)-d(1)*d(2)/9.* &
+        ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2))
+
+    dd1g12=domega/omega*dLYP12-a*b*omega*(d(2)/9.* &
+        (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt)
+
+    dd1g22=domega/omega*dLYP22-a*b*omega*(1./9.*d(2) &
+        *(1.-3.*delta-(delta-11.)*d(2)/dt)-d(1)*d(2)/9.* &
+        ((3.+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2)-2*d(1))
+
+    dd2g22=domega/omega*dLYP22-a*b*omega*(d(1)/9.* &
+        (1.-3.*delta-2*(delta-11.)*d(2)/dt)-d(1)*d(2)/9.* &
+        ((3+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2))
+
+    dd2g12=domega/omega*dLYP12-a*b*omega*(d(1)/9.* &
+        (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt)
+
+    dd2g11=domega/omega*dLYP11-a*b*omega*(1./9.*d(1) &
+        *(1.-3.*delta-(delta-11.)*d(1)/dt)-d(1)*d(2)/9.* &
+        ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2)-2*d(2))
+
+    dLYPdd(1)=-4*a/den*d(1)*d(2)/dt* &
+        (thd*dd*dt**(-fothd)/den &
+        +1./d(1)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)* &
+        (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(2)*(onzthd* &
+        d(1)**(8./3.)+d(2)**(8./3.)))+dd1g11*gam11+ &
+        dd1g12*gam12+dd1g22*gam22
+
+    dLYPdd(2)=-4*a/den*d(1)*d(2)/dt*(thd*dd*dt**(-fothd)/den &
+        +1./d(2)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)* &
+        (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(1)*(onzthd* &
+        d(2)**(8./3.)+d(1)**(8./3.)))+dd2g22*gam22+ &
+        dd2g12*gam12+dd2g11*gam11
+
+    ! second derivatives with respect to the density gradient
+    do is=1,2
+      do ix=1,3
+        dg11dd(ix,is)=2*gd(ix,is)
+        dg22dd(ix,is)=2*gd(ix,is)
+      enddo
+    enddo
+    do ix=1,3
+      dLYPgd(ix,1)=dLYP11*dg11dd(ix,1)+dLYP12*gd(ix,2)
+      dLYPgd(ix,2)=dLYP22*dg22dd(ix,2)+dLYP12*gd(ix,1)
+    enddo
+
+    EX=becke
+    EC=LYP
+    do is=1,nspin
+      dEXdd(is)=dbecdd(is)
+      dECdd(is)=dLYPdd(is)
+      do ix=1,3
+        dEXdgd(ix,is)=dbecgd(ix,is)
+        dECdgd(ix,is)=dLYPgd(ix,is)
+      enddo
+    enddo
+  end subroutine blypxc
+
+
+
+  SUBROUTINE RPBEXC( IREL, nspin, Dens, GDens, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Hammer's RPBE Generalized-Gradient-Approximation (GGA).
+    ! A revision of PBE (Perdew-Burke-Ernzerhof) 
+    ! Ref: Hammer, Hansen & Norskov, PRB 59, 7413 (1999) and
+    ! J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
+    !
+    ! Written by M.V. Fernandez-Serra. March 2004. On the PBE routine of
+    ! L.C.Balbas and J.M.Soler. December 1996.
+    ! ******** INPUT ******************************************************
+    ! INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
+    ! INTEGER nspin          : Number of spin polarizations (1 or 2)
+    ! REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
+    !                           spin electron density (if nspin=2)
+    ! REAL*8  GDens(3,nspin) : Total or spin density gradient
+    ! ******** OUTPUT *****************************************************
+    ! REAL*8  EX             : Exchange energy density
+    ! REAL*8  EC             : Correlation energy density
+    ! REAL*8  DEXDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ex)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          exchange potential
+    ! REAL*8  DECDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ec)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          correlation potential
+    ! REAL*8  DEXDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ex)/d(GradDens(i,ispin))
+    ! REAL*8  DECDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ec)/d(GradDens(i,ispin))
+    ! ********* UNITS ****************************************************
+    ! Lengths in Bohr
+    ! Densities in electrons per Bohr**3
+    ! Energies in Hartrees
+    ! Gradient vectors in cartesian coordinates
+    ! ********* ROUTINES CALLED ******************************************
+    ! EXCHNG, PW92C
+    ! ********************************************************************
+
+    INTEGER           IREL, nspin
+    real(dp)          Dens(nspin), DECDD(nspin), DECDGD(3,nspin), &
+        DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
+
+    ! Internal variables
+    INTEGER :: IS, IX
+
+    real(dp) &
+        A, BETA, D(2), DADD, DECUDD, DENMIN, &
+        DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD, &
+        DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2), &
+        DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, &
+        DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), &
+        EC, ECUNIF, EX, EXUNIF, &
+        F, F1, F2, F3, F4, FC, FX, FOUTHD, &
+        GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), &
+        H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S, &
+        T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
+
+    ! Lower bounds of density and its gradient to avoid divisions by zero
+    PARAMETER ( DENMIN = 1.D-12 )
+    PARAMETER ( GDMIN  = 1.D-12 )
+
+    ! Fix some numerical parameters
+    PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, &
+        THD=1.D0/3.D0, THRHLF=1.5D0, &
+        TWO=2.D0, TWOTHD=2.D0/3.D0 )
+
+    ! Fix some more numerical constants
+    PI = 4 * ATAN(1.D0)
+    BETA = 0.066725D0
+    GAMMA = (1 - LOG(TWO)) / PI**2
+    MU = BETA * PI**2 / 3
+    KAPPA = 0.804D0
+
+    ! Translate density and its gradient to new variables
+    IF (nspin .EQ. 1) THEN
+      D(1) = HALF * Dens(1)
+      D(2) = D(1)
+      DT = MAX( DENMIN, Dens(1) )
+      DO IX = 1,3
+        GD(IX,1) = HALF * GDens(IX,1)
+        GD(IX,2) = GD(IX,1)
+        GDT(IX) = GDens(IX,1)
+      end DO
+    ELSE
+      D(1) = Dens(1)
+      D(2) = Dens(2)
+      DT = MAX( DENMIN, Dens(1)+Dens(2) )
+      DO IX = 1,3
+        GD(IX,1) = GDens(IX,1)
+        GD(IX,2) = GDens(IX,2)
+        GDT(IX) = GDens(IX,1) + GDens(IX,2)
+      end DO
+    ENDIF
+    GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
+    GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
+    GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
+    GDMT = MAX( GDMIN, GDMT )
+
+    ! Find local correlation energy and potential
+    CALL PW92C( 2, D, ECUNIF, VCUNIF )
+
+    ! Find total correlation energy
+    RS = ( 3 / (4*PI*DT) )**THD
+    KF = (3 * PI**2 * DT)**THD
+    KS = SQRT( 4 * KF / PI )
+    ZETA = ( D(1) - D(2) ) / DT
+    ZETA = MAX( -1.D0+DENMIN, ZETA )
+    ZETA = MIN(  1.D0-DENMIN, ZETA )
+    PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
+    T = GDMT / (2 * PHI * KS * DT)
+    F1 = ECUNIF / GAMMA / PHI**3
+    F2 = EXP(-F1)
+    A = BETA / GAMMA / (F2-1)
+    F3 = T**2 + A * T**4
+    F4 = BETA/GAMMA * F3 / (1 + A*F3)
+    H = GAMMA * PHI**3 * LOG( 1 + F4 )
+    FC = ECUNIF + H
+
+    ! Find correlation energy derivatives
+    DRSDD = - (THD * RS / DT)
+    DKFDD =   THD * KF / DT
+    DKSDD = HALF * KS * DKFDD / KF
+    DZDD(1) =   1 / DT - ZETA / DT
+    DZDD(2) = - (1 / DT) - ZETA / DT
+    DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
+    DO IS = 1,2
+      DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
+      DPDD = DPDZ * DZDD(IS)
+      DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
+      DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
+      DF2DD = (- F2) * DF1DD
+      DADD = (- A) * DF2DD / (F2-1)
+      DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
+      DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
+      DHDD = 3 * H * DPDD / PHI
+      DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
+      DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
+
+      DO IX = 1,3
+        DTDGD = (T / GDMT) * GDT(IX) / GDMT
+        DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
+        DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 
+        DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
+        DFCDGD(IX,IS) = DT * DHDGD
+      end DO
+    end DO
+
+    ! Find exchange energy and potential
+    FX = 0
+    DO IS = 1,2
+      DS(IS)   = MAX( DENMIN, 2 * D(IS) )
+      GDMS = MAX( GDMIN, 2 * GDM(IS) )
+      KFS = (3 * PI**2 * DS(IS))**THD
+      S = GDMS / (2 * KFS * DS(IS))
+      !ea Hammer's RPBE (Hammer, Hansen & Norskov PRB 59 7413 (99)
+      !ea     F1 = DEXP( - MU * S**2 / KAPPA)
+      !ea     F = 1 + KAPPA * (1 - F1)
+      !ea Following is standard PBE
+      !ea     F1 = 1 + MU * S**2 / KAPPA
+      !ea     F = 1 + KAPPA - KAPPA / F1
+      !ea (If revPBE Zhang & Yang, PRL 80,890(1998),change PBE's KAPPA to 1.245)
+      F1 = DEXP( - MU * S**2 / KAPPA)
+      F = 1 + KAPPA * (1 - F1)
+
+      !       Note nspin=1 in call to exchng...
+
+      CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
+      FX = FX + DS(IS) * EXUNIF * F
+
+      !MVFS   The derivatives of F  also need to be changed for Hammer's RPBE.
+      !MVFS   DF1DD = 2 * F1 * DSDD  * ( - MU * S / KAPPA)
+      !MVFS   DF1DGD= 2 * F1 * DSDGD * ( - MU * S / KAPPA)
+      !MVFS   DFDD  = -1 * KAPPA * DF1DD
+      !MVFS   DFDGD = -1 * KAPPA * DFDGD
+
+      DKFDD = THD * KFS / DS(IS)
+      DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
+      !       DF1DD = 2 * (F1-1) * DSDD / S
+      !       DFDD = KAPPA * DF1DD / F1**2
+      DF1DD = 2* F1 * DSDD * ( - MU * S / KAPPA)
+      DFDD = -1 * KAPPA * DF1DD
+      DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
+
+      DO IX = 1,3
+        GDS = 2 * GD(IX,IS)
+        DSDGD = (S / GDMS) * GDS / GDMS
+        !         DF1DGD = 2 * MU * S * DSDGD / KAPPA
+        !         DFDGD = KAPPA * DF1DGD / F1**2
+        DF1DGD =2*F1 * DSDGD * ( - MU * S / KAPPA)
+        DFDGD = -1 * KAPPA * DF1DGD
+        DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
+      end DO
+    end DO
+    FX = HALF * FX / DT
+
+    ! Set output arguments
+    EX = FX
+    EC = FC
+    DO IS = 1,nspin
+      DEXDD(IS) = DFXDD(IS)
+      DECDD(IS) = DFCDD(IS)
+      DO IX = 1,3
+        DEXDGD(IX,IS) = DFXDGD(IX,IS)
+        DECDGD(IX,IS) = DFCDGD(IX,IS)
+      end DO
+    end DO
+  END SUBROUTINE RPBEXC
+
+
+
+  SUBROUTINE WCXC( IREL, nspin, Dens, GDens, &
+      EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+    ! *********************************************************************
+    ! Implements Wu-Cohen Generalized-Gradient-Approximation.
+    ! Ref: Z. Wu and R. E. Cohen PRB 73, 235116 (2006)
+    ! Written by Marivi Fernandez-Serra, with contributions by
+    ! Julian Gale and Alberto Garcia,
+    ! over the PBEXC subroutine of L.C.Balbas and J.M.Soler.
+    ! September, 2006.
+    ! ******** INPUT ******************************************************
+    ! INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
+    ! INTEGER nspin          : Number of spin polarizations (1 or 2)
+    ! REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
+    !                           spin electron density (if nspin=2)
+    ! REAL*8  GDens(3,nspin) : Total or spin density gradient
+    ! ******** OUTPUT *****************************************************
+    ! REAL*8  EX             : Exchange energy density
+    ! REAL*8  EC             : Correlation energy density
+    ! REAL*8  DEXDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ex)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          exchange potential
+    ! REAL*8  DECDD(nspin)   : Partial derivative
+    !                           d(DensTot*Ec)/dDens(ispin),
+    !                           where DensTot = Sum_ispin( Dens(ispin) )
+    !                          For a constant density, this is the
+    !                          correlation potential
+    ! REAL*8  DEXDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ex)/d(GradDens(i,ispin))
+    ! REAL*8  DECDGD(3,nspin): Partial derivative
+    !                           d(DensTot*Ec)/d(GradDens(i,ispin))
+    ! ********* UNITS ****************************************************
+    ! Lengths in Bohr
+    ! Densities in electrons per Bohr**3
+    ! Energies in Hartrees
+    ! Gradient vectors in cartesian coordinates
+    ! ********* ROUTINES CALLED ******************************************
+    ! EXCHNG, PW92C
+    ! ********************************************************************
+
+    INTEGER           IREL, nspin
+    real(dp)          Dens(nspin), DECDD(nspin), DECDGD(3,nspin), &
+        DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
+
+    ! Internal variables
+    INTEGER :: IS, IX
+
+    real(dp) &
+        A, BETA, D(2), DADD, DECUDD, DENMIN,  &
+        DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD, &
+        DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2), &
+        DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,  &
+        DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),  &
+        XWC, DXWCDS, CWC, &
+        EC, ECUNIF, EX, EXUNIF, &
+        F, F1, F2, F3, F4, FC, FX, FOUTHD, &
+        GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), &
+        H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S, &
+        TEN81,  &
+        T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
+
+    ! Lower bounds of density and its gradient to avoid divisions by zero
+    PARAMETER ( DENMIN = 1.D-12 )
+    PARAMETER ( GDMIN  = 1.D-12 )
+
+    ! Fix some numerical parameters
+    PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, &
+        THD=1.D0/3.D0, THRHLF=1.5D0, &
+        TWO=2.D0, TWOTHD=2.D0/3.D0 )
+    PARAMETER ( TEN81 = 10.0d0/81.0d0 )
+
+    ! Fix some more numerical constants
+    PI = 4 * ATAN(1.D0)
+    BETA = 0.066725D0
+    GAMMA = (1 - LOG(TWO)) / PI**2
+    MU = BETA * PI**2 / 3
+    KAPPA = 0.804D0
+    CWC = 0.0079325D0
+
+    ! Translate density and its gradient to new variables
+    IF (nspin .EQ. 1) THEN
+      D(1) = HALF * Dens(1)
+      D(2) = D(1)
+      DT = MAX( DENMIN, Dens(1) )
+      DO IX = 1,3
+        GD(IX,1) = HALF * GDens(IX,1)
+        GD(IX,2) = GD(IX,1)
+        GDT(IX) = GDens(IX,1)
+      end DO
+    ELSE
+      D(1) = Dens(1)
+      D(2) = Dens(2)
+      DT = MAX( DENMIN, Dens(1)+Dens(2) )
+      DO IX = 1,3
+        GD(IX,1) = GDens(IX,1)
+        GD(IX,2) = GDens(IX,2)
+        GDT(IX) = GDens(IX,1) + GDens(IX,2)
+      end DO
+    ENDIF
+    GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
+    GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
+    GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
+    GDMT = MAX( GDMIN, GDMT )
+
+    ! Find local correlation energy and potential
+    CALL PW92C( 2, D, ECUNIF, VCUNIF )
+
+    ! Find total correlation energy
+    RS = ( 3 / (4*PI*DT) )**THD
+    KF = (3 * PI**2 * DT)**THD
+    KS = SQRT( 4 * KF / PI )
+    ZETA = ( D(1) - D(2) ) / DT
+    ZETA = MAX( -1.D0+DENMIN, ZETA )
+    ZETA = MIN(  1.D0-DENMIN, ZETA )
+    PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
+    T = GDMT / (2 * PHI * KS * DT)
+    F1 = ECUNIF / GAMMA / PHI**3
+    F2 = EXP(-F1)
+    A = BETA / GAMMA / (F2-1)
+    F3 = T**2 + A * T**4
+    F4 = BETA/GAMMA * F3 / (1 + A*F3)
+    H = GAMMA * PHI**3 * LOG( 1 + F4 )
+    FC = ECUNIF + H
+
+    ! Find correlation energy derivatives
+    DRSDD = - (THD * RS / DT)
+    DKFDD =   THD * KF / DT
+    DKSDD = HALF * KS * DKFDD / KF
+    DZDD(1) =   1 / DT - ZETA / DT
+    DZDD(2) = - (1 / DT) - ZETA / DT
+    DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
+    DO IS = 1,2
+      DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
+      DPDD = DPDZ * DZDD(IS)
+      DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
+      DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
+      DF2DD = (- F2) * DF1DD
+      DADD = (- A) * DF2DD / (F2-1)
+      DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
+      DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
+      DHDD = 3 * H * DPDD / PHI
+      DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
+      DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
+
+      DO IX = 1,3
+        DTDGD = (T / GDMT) * GDT(IX) / GDMT
+        DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
+        DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 
+        DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
+        DFCDGD(IX,IS) = DT * DHDGD
+      end DO
+    end DO
+
+    ! Find exchange energy and potential
+    FX = 0
+    DO IS = 1,2
+      DS(IS)   = MAX( DENMIN, 2 * D(IS) )
+      GDMS = MAX( GDMIN, 2 * GDM(IS) )
+      KFS = (3 * PI**2 * DS(IS))**THD
+      S = GDMS / (2 * KFS * DS(IS))
+      !
+      ! For PBE: 
+      !
+      !       x = MU * S**2
+      !       dxds = 2*MU*S
+      !
+      ! Wu-Cohen form:
+      !
+      XWC= TEN81 * s**2 + (MU- TEN81) * &
+          S**2 * exp(-S**2) + log(1+ CWC * S**4)
+      DXWCDS = 2 * TEN81 * S + (MU - TEN81) * exp(-S**2) * &
+          2*S * (1 - S*S) + 4 * CWC * S**3 / (1 + CWC * S**4) 
+
+      F1 = 1 +  XWC / KAPPA
+      F = 1 + KAPPA - KAPPA / F1
+      !
+      !       Note nspin=1 in call to exchng...
+      !
+      CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
+      FX = FX + DS(IS) * EXUNIF * F
+
+      DKFDD = THD * KFS / DS(IS)
+      DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
+      DF1DD = DXWCDS * DSDD / KAPPA 
+      DFDD = KAPPA * DF1DD / F1**2
+      DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
+
+      DO IX = 1,3
+        GDS = 2 * GD(IX,IS)
+        DSDGD = (S / GDMS) * GDS / GDMS
+        DF1DGD = DXWCDS * DSDGD / KAPPA
+        DFDGD = KAPPA * DF1DGD / F1**2
+        DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
+      end DO
+    end DO
+    FX = HALF * FX / DT
+
+    ! Set output arguments
+    EX = FX
+    EC = FC
+    DO IS = 1,nspin
+      DEXDD(IS) = DFXDD(IS)
+      DECDD(IS) = DFCDD(IS)
+      DO IX = 1,3
+        DEXDGD(IX,IS) = DFXDGD(IX,IS)
+        DECDGD(IX,IS) = DFCDGD(IX,IS)
+      end DO
+    end DO
+