subroutine exchen(d,s,ex)
c  gga91 exchange energy for a spin-unpolarized electronic system
c  input d   : density
c  input s   : abs(grad d)/(2*kf*d)
c  output ex : exchange energy per electron 
      implicit double precision (a-h,o-z)
      data a1,a2,a3,a4/0.19645d0,0.27430d0,0.15084d0,100.d0/
      data ax,a,b1/-0.7385588d0,7.7956d0,0.004d0/
      data thrd/0.333333333333d0/
      fac = ax*d**thrd
      s2 = s*s
      s4 = s2*s2
      p0 = 1.d0/dsqrt(1.d0+a*a*s2)
      p1 = dlog(a*s+1.d0/p0)
      p2 = dexp(-a4*s2)
      p3 = 1.d0/(1.d0+a1*s*p1+b1*s4)
      p4 = 1.d0+a1*s*p1+(a2-a3*p2)*s2
      f = p3*p4
      ex = fac*f
      return
      end