subroutine LSD2D(rs,z,ec,vcup,vcdown)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc computes the correlation energy of the 2D electron gas according
cc to Eqs. (3), (4) and Table II of Attaccalite et al., 
cc Phys. Rev. Lett. 88, 256601 (2002); 91, 109902(E) (2003), and
cc the corresponding LSD correlation potential according
cc to Eqs. (17)-(23) and Table I of Gori-Giorgi et al., 
cc Int. J. Quantum Chem. 91, 126 (2003) [cond-mat/0110444]
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc energies in Hartree atomic units
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc ec -> correlation energy
cc vcup -> LSD correlation potential for spin-up electrons
cc vcdown -> LSD correlation potential for spin-down electrons
cc (rs,z) -> (r_s,\zeta)
      implicit none
      real*8 rs,z,ec,vcup,vcdown
      real*8 ex6,ecdrs,ecdz,cF,cFd,ax,z4
      real*8 pi,A,B,C,D,E,F,G,H
      real*8 alpha0,alpha0d,alpha1,alpha1d,alpha2,alpha2d,bet
      real*8 am1,am32,am2_0,am2_1,am2_2

      pi=acos(-1.d0)
ccc alpha_0(rs)
      A=-0.1925d0
      B=0.0863136d0
      C=0.0572384d0
      E=1.0022d0
      F=-0.02069d0
      G=0.33997d0
      H=1.747d-2
      D=-A*H
      alpha0=A+(B*rs+C*rs**2+D*rs**3)*log(1+1/(E*rs+
     $     F*rs**(3.d0/2.d0)+G*rs**2+H*rs**3))
ccc d/drs alpha_0(rs)
      alpha0d=(B+2*C*rs+3*D*rs**2)*log(1+1/(E*rs+F*sqrt(rs)**3+G*rs**2
     $     +H*rs**3))-(B*rs+C*rs**2+D*rs**3)/(E*rs+F*sqrt(rs)**3+G*rs**2
     $     +H*rs**3)**2*(E+3.D0/2.D0*F*sqrt(rs)+2*G*rs+3*H*rs**2)/(1+1/
     $     (E*rs+F*sqrt(rs)**3+G*rs**2+H*rs**3))

ccc coefficients for the large rs expansion
      am1=C/H+A*G/H             !1/rs
      am32=A*F/H                !1/rs^(3/2)
      am2_0=-1./H**2*(-A*E*H+A*G**2-B*H+C*G) !1/rs^2

ccc alpha_1(rs)
      A=0.117331d0
      B=-3.394d-2
      C=-7.66765d-3
      E=0.4133d0
      F=0.d0
      G=6.68467d-2
      H=7.799d-4
      D=-A*H      
      alpha1=A+(B*rs+C*rs**2+D*rs**3)*log(1.+1./(E*rs+
     $     F*rs**(3.d0/2.d0)+G*rs**2+H*rs**3))
ccc d/drs alpha_1(rs)
      alpha1d=(B+2*C*rs+3*D*rs**2)*log(1+1/(E*rs+F*sqrt(rs)**3+G*rs**2
     $     +H*rs**3))-(B*rs+C*rs**2+D*rs**3)/(E*rs+F*sqrt(rs)**3+G*rs**2
     $     +H*rs**3)**2*(E+3.D0/2.D0*F*sqrt(rs)+2*G*rs+3*H*rs**2)/(1+1/
     $     (E*rs+F*sqrt(rs)**3+G*rs**2+H*rs**3))

ccc large rs epansion (1/rs cancels with ex)
      am2_1=-1.d0/H**2*(-A*E*H+A*G**2-B*H+C*G) !1/rs^2

ccc alpha_2(rs)
      A=0.0234188d0
      B=-0.037093d0
      C=0.0163618d0
      E=1.424301d0
      F=0.d0
      G=0.d0
      H=1.163099d0
      D=-A*H
      alpha2=A+(B*rs+C*rs**2+D*rs**3)*log(1+1/(E*rs+
     $     F*rs**(3.d0/2.d0)+G*rs**2+H*rs**3))
ccc d/drs alpha_2(rs)
      alpha2d=(B+2*C*rs+3*D*rs**2)*log(1+1/(E*rs+F*sqrt(rs)**3+G*rs**2
     $     +H*rs**3))-(B*rs+C*rs**2+D*rs**3)/(E*rs+F*sqrt(rs)**3+G*rs**2
     $     +H*rs**3)**2*(E+3.D0/2.D0*F*sqrt(rs)+2*G*rs+3*H*rs**2)/(1+1/
     $     (E*rs+F*sqrt(rs)**3+G*rs**2+H*rs**3))

ccc large rs epansion (1/rs cancels with ex)
      am2_2=-1.d0/H**2*(-A*E*H+A*G**2-B*H+C*G) !1/rs^2

      bet=1.3386d0                !beta
      z4=z**4
      ax=4./(3.*pi*sqrt(2.))    !ax
      cF=(1.d0+z)**(3.d0/2.d0)+(1.d0-z)**(3.d0/2.d0)-(2.d0+3.d0/4.d0*z*z
     $     +3.d0/64.d0*z4)                  !\cal{F}
      ex6=-cF*ax/rs

ccc d/dz cF
      cFd=3.d0/2.d0*(sqrt(1.d0+z)-sqrt(1.d0-z))-3.d0/2.d0*z-
     $     3.d0/16.d0*z**3

ccc correlation energy and its derivatives w.r.t rs and z ccccccccccccc
ccc switch to asymptotic expansion when rs > 4000
      if(rs.lt.4.d3) then
         ec=(exp(-bet*rs)-1.d0)*ex6+alpha0+alpha1*z*z+alpha2*z4
ccc d/drs ec
         ecdrs=ax*cF/rs/rs*(exp(-bet*rs)*(1.d0+bet*rs)-1.d0)+alpha0d+
     $        alpha1d*z*z+alpha2d*z4
ccc d/dz ec 
         ecdz=ax/rs*(1.d0-exp(-bet*rs))*cFd+2.d0*alpha1*z+4.*
     $        alpha2*z*z*z
      else
ccc ec for rs > 4000 (asymptotic exp.)
         ec=ax/rs*((1.d0+z)**(3.d0/2.d0)+(1.-z)**(3.d0/2.d0))-2.d0*ax/rs
     $        +am1/rs+am32/rs**(3.d0/2.d0)+(am2_0+z**2*am2_1+z**4*am2_2)
     $        /rs**2
ccc d/drs ec for rs > 4000
         ecdrs=-(ax*((1.d0+z)**(3.d0/2.d0)+(1.d0-z)**(3.d0/2.d0))-2.d0
     $        *ax+am1)/rs**2-3.d0/2.d0*am32/rs**(5.d0/2.d0)-2.d0/rs**3
     $        *(am2_0+z**2*am2_1+z**4*am2_2)
ccc d/dz ec for rs > 4000
         ecdz=ax*3.d0/2.d0*((1.d0+z)**(1.d0/2.d0)-(1.d0-z)**(1.d0/2.d0)
     $        )/rs+(2.d0*z*am2_1+4.d0*z**3*am2_2)/rs**2
      endif
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc corr. pot. for spin-up electrons
      vcup=ec-rs/2.d0*ecdrs-(z-1.d0)*ecdz
ccccccccccccccccccccccccccccccccccccc
ccc corr. pot. for spin-down electrons
      vcdown=ec-rs/2.d0*ecdrs-(z+1.d0)*ecdz
ccccccccccccccccccccccccccccccccccccc
      return
      end