1 files changed, 634 insertions, 0 deletions

diff --git a/src/labat/resten.f b/src/labat/resten.f
new file mode 100644
index 0000000..8314dad
--- /dev/null
+++ b/src/labat/resten.f
@@ -0,0 +1,634 @@
+c *********************************************************************
+      subroutine scheq(zz,en,lambda,nofl,kkk,mess,scf,thresh,iflag,npr
+     1,de)
+      implicit double precision (a-h,o-z)
+c     subroutine scheq
+c     compute energy eigenvalue and wave function
+c      originally written by sherwood skillman
+c      rca laboratories, princeton, new jersey, spring 1961
+c      modified by frank herman, summer 1961
+c      further modified by richard kortum and paul kelly, lockheed
+c       research laboratories, palo alto, california,  summer 1962
+c further modified by r.m. nieminen, otaniemi,spring 1978
+c and last but not least modified by m.p. winter 1979
+      dimension qq(561),p(5),q(5),t(5),d(5)
+      common/sc/gr(561),r(561),snlo(561),nbl
+      common/pot/v(561)
+      imp=1
+      nb1=nbl+1
+      nmb=nb1/5
+      if (de.lt.1.d-10)imp=0
+      if(imp.eq.1)de=dmax1(de,-e/500.d0)
+      eg=0.d0
+       z=zz
+       lam=lambda
+       nn=nofl
+       mesh=mess
+       c=scf
+          i0=0
+      many = 150
+      iflag=0
+   73 e=en
+      morev=0
+      lessv=0
+      emore=0.0d0
+      eless=0.0d0
+      more=0
+      less=0
+      mp=0
+      nprint=0
+      lamm=lam-1
+      lamp=lam+1
+      xlp=lamp
+      ndcr=nn-lamp
+      b=lam*lamp
+       oc=r(2)
+      h=oc
+      hsq=h*h
+      b3=(v(3)-v(2))/h-z/hsq
+      y=h+h
+      flps=4*lam+6
+      slpt=6*lam+12
+      elpt=8*lam+20
+      a1=-z/xlp
+      ysq=y*y
+       b1=-z-z
+      ab1=a1*b1
+      ab3=a1*b3
+c     raise h and y to lam+1
+      htl=h
+      ytl=y
+      if(lam)77,1102,1100
+ 1100 do 1101 i=1,lam
+      htl=htl*h
+ 1101  ytl=ytl*y
+ 1102 h1=hsq
+      bohs=b/hsq
+      boh=b1/h
+      bth=b3*h
+      bq3=bohs+boh+bth
+      bq4=bohs/4.d0+boh/2.d0+bth+bth
+      epl=8+lam
+      fpl=5+lam
+      xifc=c*.21701389d-4
+c     start outward integration
+   10 nprint=nprint+1
+      eps =e-eg
+      eg =e
+      if(many-nprint)  900,75,75
+900   if(mp.eq.2)go to 300
+      write (6,901)nn,lam ,z
+      call flush(6)
+  901 format (21h    no convergence on,i4,i1,f4.0)
+      iflag=1
+      en=0.d0
+      npr=nprint
+      return
+  77  nstop=7
+ 777  write(66,770) nstop
+ 770  format(' error in scheq',i6)
+      en=0.d0
+      iflag=1
+      npr=nprint
+      return
+   75 do 11 i=1,mesh
+   11 snlo(i)=0.0d0
+      if(nprint-1)  77,500,529
+  500 continue
+      do 502  i= 4,mesh
+      qq(i) = v(i)+b/(r(i)*r(i))-e
+  502 continue
+  516 m= mesh
+      do 520 i=4,mesh
+      if(qq(m)) 519,520,520
+  519 ik=m+1
+      go to 525
+  520 m=m-1
+      write(66,562)lam,qq(m+1),e
+562   format(/'  l: ',i5,'  q:  ',g15.6,'  e:  ',g15.6/)
+  521 nstop =521
+c     q is everywhere positive
+      go to 777
+  525 if(mesh-ik) 526,526,535
+  526 eps = qq(mesh-nbl)
+      e = e+eps
+  529 continue
+      do 530 i=4,mesh
+  530 qq(i) = qq(i)-eps
+      go to 516
+  535 continue
+   14 ncross=0
+      sig=1.0d0
+      h=oc
+      y=h+h
+c     b= lam*(lam+1)
+c     b1= -2.d0*z
+      b2=3.0d0*z/h-e+2.0d0*v(2)-v(3)
+c     b3=(v(3)-v(2))/h -z/hsq
+c     a1= -z/(lam+1)
+      a2=(ab1+b2)/flps
+c     a2=(a1*b1+b2)/(4*lam+6)
+      a3=(a2*b1+a1*b2+b3)/slpt
+c     a3=(a2*b1+a1*b2+b3)/(6*lam+12)
+      a4=(a3*b1+a2*b2+ab3)/elpt
+c     a4=(a3*b1+a2*b2+a1*b3)/(8*lam+20)
+      p(3)=(1.0d0+h*(a1+h*(a2+h*(a3+h*a4))))*htl
+c     p(3)=(1.0d0+a1*h+a2*h**2+a3*h**3+a4*h**4)*h**(xl+1.0d0)
+      p(4)=(1.0d0+y*(a1+y*(a2+y*(a3+y*a4))))*ytl
+c     p(4)=(1.0d0+a1*y+a2*y**2+a3*y**3+a4*y**4)*y**(xl+1.0d0)
+      q(3)=bq3+b2
+c     q(3)=(b+b1*h+b2*h**2+b3*h**3)/h**2
+      q(4)=bq4+b2
+c     q(4)=(b+b1*y+b2*y**2+b3*y**3)/y**2
+      snlo(2)=p(3)
+      snlo(3)=p(4)
+      i=3
+       dx=oc
+      h1=h**2
+      h2=h1/12.0d0
+      t(3)=p(3)*(1.0d0-h2*q(3))
+      t(4)=p(4)*(1.0d0-h2*q(4))
+      d(4)=t(4)-t(3)
+      ncount=3
+      nint=2
+   15 i=i+1
+c     if end of mesh is reached, modify trial eigenvalue
+      if(i-mesh) 16,3,3
+    3 if(ndcr-ncross) 32,33,33
+c     return to begininning of outward integration if necessary
+   16 q(5) =qq(i)
+      if(ik-i)  29,29,21
+   21 d(5)=d(4)+h1*q(4)*p(4)
+      t(5)=d(5)+t(4)
+      if(1.0d0-abs (h2*q(5))) 3,3,501
+  501 p(5)=t(5)/(1.0d0-h2*q(5))
+      snlo(i) = p(5)
+       if(sig) 211,77,212
+  211  if(p(5))  23,23,22
+  212  if(p(5)) 22,23,23
+   22 ncross=ncross+1
+c     count changes in sign
+      sig=-sig
+   23 ncount=ncount+1
+      if(7-ncount)77,24,25
+   24 ncount=2
+   25 nint=nint+1
+      if(nbl-nint)77,26,27
+   26 dx=dx+dx
+      h=dx
+      h1=h**2
+      h2=h1/12.0d0
+      nint=0
+      t(5)=p(5)*(1.0d0-h2*q(5))
+      t(3)=p(3)*(1.0d0-h2*q(3))
+      d(5)=t(5)-t(3)
+   27 do 28 k=1,4
+      p(k)=p(k+1)
+      t(k)=t(k+1)
+      d(k)=d(k+1)
+   28 q(k)=q(k+1)
+      go to 15
+   29 if(ncount-2)77,30,21
+   30 if(nint-4)21,21,31
+c     matching radius has been reached going out
+c     if ndcr not equal to ncross,modify trial eigenvalue
+   31 eigen=e
+      if(ndcr-ncross) 32,35,33
+   32 more=1
+      mp=0
+c     too many crossings, increase absf(e)
+      morev=morev+1
+       if(morev-1) 50,53,52
+   50 nstop=50
+      go to 777
+   52 if(e -emore) 53,54,54
+   53 emore=e
+   54  if (less) 55,56,64
+   55 nstop=55
+      go to 777
+   56 e=1.25d0*eg
+      go to 10
+   33 less=1
+      mp=0
+c     too few  crossings, decrease abs (e)
+      lessv=lessv+1
+       if(lessv-1) 57,60,59
+   57 nstop=57
+      go to 777
+   59 if(eless- e) 60,61,61
+   60 eless=e
+   61 if(more) 62,63,64
+   62 nstop=62
+      go to 777
+   63 e=0.75d0*eg
+      go to 10
+   64 e=0.5d0*(emore+eless)
+      go to 10
+   35  if(abs (snlo(i-1))-abs (snlo(i-2))) 351,354,354
+c     check to see that wave is in the damped region (absolute value
+c     decreasing and signs alike)
+  351 if (p(5)) 352,21 ,353
+  352 if(snlo(i-2)) 401,21,21
+  353 if(snlo(i-2)) 21,21,401
+  354 if(1.0d+25 -abs (p(5)))33,33,21
+c     large absolute value of p in what should be the damped region
+c       indicates too few peaks, decrease absf(e)
+c     now ndcr = ncross and matching radius lies in damped region
+  401 imatch=i-2
+      xmatch=r(i-2)
+      ppout=(t(4)-t(2)-0.5d0*(p(4)-p(2)))/h
+      s2=ppout/p(3)
+c     integration is by 8 applications of newton-cotes closed
+c       quadrature for five intervals on each block
+c       xifc =(5*h(block 1)/288)/2 ,h(1) =0.0025*scale factor
+      sum1=0.0d0
+      xif=xifc
+      i=1
+      value=0.0d0
+   36 mm=nmb
+      sum2=0.0d0
+      xif=xif+xif
+   37 y=value
+      value=snlo(i+5)**2
+      sum2=sum2+19.d0*(value+y)+75.d0*(snlo(i+4)**2+snlo(i+1)**2)
+     1 +50.0d0*(snlo(i+2)**2+snlo(i+3)**2)
+      i=i+5
+       if (imatch-i)  77,39,371
+  371  mm=mm-1
+      if(mm)77,38,37
+   38 sum1=sum2*xif+sum1
+       go to 36
+   39 sum1= sum1+sum2*xif
+   40 s1=sum1/p(3)**2
+      pmatch=p(3)
+      if(nn-1)77,41,42
+   41 xinw=epl* xmatch
+c     for n =1, start inward integration at(8+lam)*xmatch or x max
+      go to 421
+   42 xinw=fpl*xmatch
+c     for n not=1, start at (5+lam)*xmatch or x max (end of mesh)
+  421 do 44  i= nb1,mesh,nbl
+      if(xinw-r(i))  43,43,44
+   43 kkk =i
+      go to 45
+   44 continue
+      kkk =mesh
+   45 i =kkk
+      dx =r(i-1)-r(i)
+      h =dx
+      xif=0.17361111d-1*dx
+      hsq=h*h
+      hsq12=hsq/12.0d0
+      q(3)= qq(i)
+      p(3)= exp (-r(i)*sqrt (q(3)))
+  402 sum3=p(3)/q(3)
+      i=i-1
+      q(4) =qq(i)
+  404 p(4)=exp (-r(i)*sqrt (q(4)))
+      if (abs (p(4))-1.0d-35) 4041,4041,405
+ 4041 kkk=kkk -nbl
+      if(kkk -imatch) 4042,4042,45
+ 4042 write (5,4043)z ,nn,lam,kkk
+ 4043 format  (6hoat z=,f6.0,6h  nl =,i3,i1,7h  kkk =,i5,22h is less tha
+     1n imatch =,i5,44h inward integration will be tried at kkk+nbl)
+      kkk =kkk +nbl
+      p(4) = 1.5d-35
+      p(3) = 1.0d-35
+  405 if(pmatch)102,77,103
+  102 p(3)=-p(3)
+      p(4)=-p(4)
+  103 snlo(i+1)=p(3)
+      snlo(i)=p(4)
+      t(3)=p(3)*(1.0d0-hsq12 *q(3))
+      t(4)=p(4)*(1.0d0-hsq12*q(4))
+      d(4)=t(4)-t(3)
+  104 do 106 m=2,nbl
+      i=i-1
+      q(5) =qq(i)
+      d(5)=hsq*q(4)*p(4)+d(4)
+      t(5)=d(5)+t(4)
+      p(5)=t(5)/(1.0d0-hsq12*q(5))
+      if(i-imatch+1)77,200,105
+  105 snlo(i)=p(5)
+      do 106 k=1,4
+      p(k)=p(k+1)
+      t(k)=t(k+1)
+      d(k)=d(k+1)
+  106 q(k)=q(k+1)
+      q(5) =qq(i-2)
+      d(5)=hsq*q(4)*p(4)+d(4)
+      t(5)=d(5)+t(4)
+      p(5)=t(5)/(1.0d0-hsq12*q(5))
+      p(5)=1.09375d0*p(4)+0.2734375d0*p(5)-0.546875d0*p(3)+0.21875d0*
+     1 p(2)-0.0390625d0*p(1)
+      i=i-1
+      dx=dx/2.0d0
+      q(5) =qq(i)
+      h=dx
+      hsq=h*h
+      hsq12=hsq/12.0d0
+      t(5)=p(5)*(1.0d0-hsq12*q(5))
+      t(4)=p(4)*(1.0d0-hsq12*q(4))
+      d(5)=t(5)-t(4)
+      snlo(i)=p(5)
+      do 107 l=1,4
+      p(l)=p(l+1)
+      t(l)=t(l+1)
+      d(l)=d(l+1)
+  107 q(l)=q(l+1)
+      go to 104
+c     matching radius has been reached coming in
+  200 k=kkk
+      value=snlo(k)**2
+      go to 202
+ 2001  continue
+  201 sum3=sum3+xif*sum4
+      xif =xif*0.5d0
+  202 mm=nmb
+      sum4 =0.0d0
+  203 y=value
+      value=snlo(k-5)**2
+      sum4=sum4+19.d0*(value+y)+75.d0*(snlo(k-1)**2+snlo(k-4)**2)
+     1+50.0d0*(snlo(k-2)**2+snlo(k-3)**2)
+      k=k-5
+       if(k-imatch) 77,2031,2030
+ 2030  mm=mm-1
+       if(mm)  77,2001,203
+ 2031 sum3=sum3+xif*sum4
+  204 s3=sum3/p(4)**2
+      ppin=(t(5)-t(3)-0.5d0*(p(5)-p(3)))/h
+      s4=ppin/p(4)
+      fe=s4-s2
+      if(abs(fe).lt.1.d-08)go to 300
+c         write(66,4711)mp,nprint,e,fe
+          ml=more*less
+4711  format(' mp',i2,' iter: ',i4,' e: ',e12.6,' fe: ',e12.6)
+      if(mp.eq.2)go to 9120
+      if(mp.eq.1)go to 9100
+      mp=1
+      if(imp.eq.0)de=-e/10.d0
+          de=dmin1(de,-e/10.d0)
+          if(i0.gt.0)de=de/5**i0
+          i0=i0+1
+        if(ml.gt.0)de=abs(emore-eless)/10.d0
+      if(e.gt.-1.d-35)de=1.d-7
+      eold=e
+      feold=fe
+      if(fe.gt.0)e=e-de
+      if(fe.lt.0)e=e+de
+      go to 10
+9100  if(fe*feold.lt.0)go to 9110
+      if(fe.gt.0.d0.and.imp.eq.0) de=-e/10.d0
+          if(i0.gt.1.and.ml.gt.0)de=abs(emore-eless)/10.d0
+      eold=e
+      feold=fe
+9101  if(fe.lt.0d0)e=e+de
+      if(fe.gt.0d0)e=e-de
+      go to 10
+9110  continue
+9111  mp=2
+      fe1=feold
+      fe2=fe
+      e1=eold
+      e2=e
+      eold=e
+      go to 9130
+9120  if(fe1*fe.lt.0d0)e2=e
+      if(fe1*fe.lt.0d0)fe2=fe
+      if(fe2*fe.lt.0d0)e1=e
+      if(fe2*fe.lt.0d0)fe1=fe
+      eold=e
+6174  format(' e1: ',e12.6,' fe1: ',e12.6,' e2: ',e12.6,
+     j' fe2: ',e12.6)
+9130  e=(e1*fe2-e2*fe1)/(fe2-fe1)
+      if(abs((e2-e1)/e).lt.thresh)go to 300
+      if(abs((e-eold)/e).lt.1d-06)go  to 300
+      go to 10
+  300  pop=pmatch/p(4)
+      do 302 j=imatch,kkk
+  302 snlo(j)=snlo(j)*pop
+      sum1=0.0d0
+      j=1
+      xif=xifc
+      value=0.0d0
+  303 mm=nmb
+      xif=xif+xif
+      sum2=0.0d0
+  304 y=value
+      value=snlo(j+5)**2
+      sum2=sum2+19.d0*(value+y)+75.d0*(snlo(j+4)**2+snlo(j+1)**2)
+     1+50.d0*(snlo(j+2)**2+snlo(j+3)**2)
+      j=j+5
+      mm=mm-1
+      if(mm)77,305,304
+  305 sum1=sum1+xif*sum2
+      if(kkk-j)77,307,303
+  307 c1=sqrt (sum1)
+      if(snlo(3))308,77,310
+  308 c1=-c1
+c
+c include tail
+c
+ 310  tail=snlo(mesh)**2/(2.d0*sqrt(abs(e)))
+      c1=c1+tail
+      do 311 i=1,kkk
+  311 snlo(i)=snlo(i)/c1
+      en=e
+      if(abs(en).lt.1.d-6) iflag=1
+      npr=nprint
+      return
+      end
+c **********************************************************************
+      subroutine simpsh(ff,s)
+      implicit double precision (a-h,o-z)
+c integrates the matrix ff on the herman-skillman mesh
+      common/mess/wj(301),dx,nblock,jblock
+      common/sc/gr(561),r(561),snlo(561),nbl
+      dimension ff(561)
+      nb1=nbl+1
+      s=0.d0
+      deltax=dx
+      do 10 j=1,jblock
+      do 11 jk=1,nb1
+      i=(j-1)*nbl+jk
+  11  s=s+ff(i)*wj(jk)*deltax
+  10  deltax=2.d0*deltax
+      return
+      end
+c **********************************************************************
+      subroutine scrhs(roo,u,z,q,n)
+      implicit double precision (a-h,o-z)
+c evaluates the screened poisson integral
+c in the herman-skillman mesh
+      dimension roo(561),u(561),ss(561),apu(561)
+      common/sc/gr(561),r(561),snlo(561),nbl
+      pi=4.d0*atan(1.d0)
+      c=4.d0*pi
+      do 11 i=1,n
+  11  ss(i)=-c*roo(i)*r(i)-q**2*u(i)
+      do 20 i=1,n
+      do 30 j=1,n
+      aa=gr(i)/gr(j)
+      if(j.gt.i) aa=1.d0/aa
+      ab=gr(i)*gr(j)
+ 30   apu(j)=(aa-ab)*ss(j)
+      call simpsh(apu,s)
+ 20   u(i)=-s/(2.d0*q)-z*gr(i)
+      return
+      end
+c **********************************************************************
+      function exc(r,s,ixc)
+      implicit double precision (a-h,o-z)
+c          ceperley alder perdew zunger
+c          or von barth hedin
+      data gp,bp1,bp2,gf,bf1,bf2/-.1423d0,1.0529d0,.3334d0,-.0843d0,
+     j1.3981d0,.2611d0/
+      data ap,bp,cp,dp,af,bf,cf,df/.0311d0,-.048d0,.002d0,-.0116d0,
+     j.01555d0,-.0269d0,.0007d0,-.0048d0/
+          if(r.lt.1.d-25)go to 100
+      pi=4.d0*atan(1.d0)
+       if(s.gt.0.99999999d0)s=0.99999999d0
+       if(s.lt.-0.99999999d0)s=-0.99999999d0
+          sinv=(4.d0*pi*r/3.d0)**(1.d0/3.d0)
+          rs=1.d0/sinv
+          if(ixc.eq.1)go to 200
+      if(rs.lt.1.d0)go to 10
+      srs=sqrt(rs)
+      excp=gp/(1.d0+bp1*srs+bp2*rs)
+      excf=gf/(1.d0+bf1*srs+bf2*rs)
+      go to 20
+10    aa=log(rs)
+      excp=ap*aa+bp+cp*rs*aa+dp*rs
+      excf=af*aa+bf+cf*rs*aa+df*rs
+20    f=((1.d0+s)**(4.d0/3.d0)+(1.d0-s)**(4.d0/3.d0)-2.d0)/
+     1  (2.d0**(4.d0/3.d0)-2.d0)
+      excf=excf-.5772521d0/rs
+      excp=excp-.4581653d0/rs
+      exc=excp+f*(excf-excp)
+      return
+100   exc=0.d0
+      return
+200   x=s/2.d0+0.5d0
+      d43=4.d0/3.d0
+      rsf=rs/75.d0
+      rsf2=rsf*rsf
+      rsf3=rsf2*rsf
+      rsp=rs/30.d0
+      rsp2=rsp*rsp
+      rsp3=rsp2*rsp
+      fcf=(1.d0+rsf3)*log(1.d0+1.d0/rsf)+0.5d0*rsf-rsf2-1.d0/3.d0
+      fcp=(1.d0+rsp3)*log(1.d0+1.d0/rsp)+0.5d0*rsp-rsp2-1.d0/3.d0
+      epscp=-.0504d0*fcp
+      epscf=-.0254d0*fcf
+      epsxp=-.91633059d0/rs
+      cny=5.1297628d0*(epscf-epscp)
+      aa=.5d0**(1.d0/3.d0)
+      if(x.lt..000001d0) x=.000001d0
+      if(x.gt..999999d0) x=.999999d0
+      fx=(x**d43+(1.d0-x)**d43-aa)/(1.d0-aa)
+      exc=epsxp+epscp+fx*(cny+4.d0/3.d0*epsxp)/5.1297628d0
+      exc=exc/2.d0
+      return
+      end
+c **********************************************************************
+      function uxc(r,s,ispin,ixc)
+      implicit double precision (a-h,o-z)
+c          ceperley alder perdew zunger
+      data gp,bp1,bp2,gf,bf1,bf2/-.1423d0,1.0529d0,.3334d0,-.0843d0,
+     j1.3981d0,.2611d0/
+      data ap,bp,cp,dp,af,bf,cf,df/.0311d0,-.048d0,.002d0,-.0116d0,
+     j.01555d0,-.0269d0,.0007d0,-.0048d0/
+          if(r.lt.1.d-35)go to 100
+      pi=4.d0*atan(1.d0)
+          sinv=(4.d0*pi*r/3.d0)**(1.d0/3.d0)
+          rs=1.d0/sinv
+          if(ixc.eq.1)go to 200
+      if(rs.lt.1.d0)go to 10
+      srs=sqrt(rs)
+      excp=gp/(1.d0+bp1*srs+bp2*rs)
+      excf=gf/(1.d0+bf1*srs+bf2*rs)
+      uxcp=excp*(1.d0+7.d0/6.d0*bp1*srs+4.d0/3.d0*bp2*rs)
+     j/(1.d0+bp1*srs+bp2*rs)
+      uxcf=excf*(1.d0+7.d0/6.d0*bf1*srs+4.d0/3.d0*bf2*rs)
+     j/(1.d0+bf1*srs+bf2*rs)
+      go to 20
+10    aa=log(rs)
+      excp=ap*aa+bp+cp*rs*aa+dp*rs
+      excf=af*aa+bf+cf*rs*aa+df*rs
+      uxcp=ap*aa+(bp-ap/3.d0)+2.d0/3.d0*cp*rs*aa+(2.d0*dp-cp)*rs/3.d0
+      uxcf=af*aa+(bf-af/3.d0)+2.d0/3.d0*cf*rs*aa+(2.d0*df-cf)*rs/3.d0
+20    f=((1.d0+s)**(4.d0/3.d0)+(1.d0-s)**(4.d0/3.d0)-2.d0)/
+     j (2.d0**(4.d0/3.d0)-2.d0)
+      ddf=4.d0/3.d0*((1.d0+s)**(1.d0/3.d0)-(1.d0-s)**(1.d0/3.d0))
+     j /(2.d0**(4.d0/3.d0)-2.d0)
+      uxc=uxcp+f*(uxcf-uxcp)+(excf-excp)*(3.d0-2.d0*ispin-s)*ddf
+     j -.6108871d0/rs*(1.d0+(3.d0-2.d0*ispin)*s)**(1.d0/3.d0)
+      return
+100   uxc=0.d0
+      return
+200   x=s/2.d0+0.5d0
+      d43=4.d0/3.d0
+      xx=0.5d0-s/2.d0
+      rsf=rs/75.d0
+      rsf2=rsf*rsf
+      rsf3=rsf2*rsf
+      rsp=rs/30.d0
+      rsp2=rsp*rsp
+      rsp3=rsp2*rsp
+      fcf=(1.d0+rsf3)*log(1.d0+1.d0/rsf)+0.5d0*rsf-rsf2-1.d0/3.d0
+      fcp=(1.d0+rsp3)*log(1.d0+1.d0/rsp)+0.5d0*rsp-rsp2-1.d0/3.d0
+      epscp=-.0504d0*fcp
+      epscf=-.0254d0*fcf
+      epsxp=-.91633059d0/rs
+      cny=5.1297628d0*(epscf-epscp)
+      aa=.5d0**(1.d0/3.d0)
+      if(x.lt..000001d0) x=.000001d0
+      if(xx.lt..000001d0) xx=.000001d0
+      if(x.gt..999999d0) x=.999999d0
+      if(xx.gt..999999d0) xx=.999999d0
+      ars=-1.22177412d0/rs+cny
+      brs=-0.0504d0*log(1.d0+30.d0/rs)-cny
+      trx1=(2.d0*x)**(1.d0/3.d0)
+      trx2=(2.d0*xx)**(1.d0/3.d0)
+      if(ispin.eq.1)vxc=ars*trx1+brs
+      if(ispin.eq.2)vxc=ars*trx2+brs
+      uxc=vxc/2.d0
+      return
+      end
+c ************************************************************
+      subroutine schrhs(v,k,l,ukl)
+      implicit double precision (a-h,o-z)
+c integrates the radial schrodinger equation
+c on the herman-skillman mesh
+      double precision k
+      dimension v(561),ukl(561),y(561)
+      common/mess/wj(301),dx,nblock,jblock
+      common/sc/gr(561),r(561),snlo(561),nbl
+      sd(i)=(2.d0*v(i)-k**2+l*(l+1)/r(i)**2)*ukl(i)
+      sdd(i)=(2.d0*v(i)-k**2+l*(l+1)/r(i)**2)*y(i)
+c initial values have been stored at ukl(1)...ukl(4)
+       nb1=nbl+1
+      hh=dx
+      do 10 i=5,nb1
+      ukl(i)=ukl(i-1)+ukl(i-3)-ukl(i-4)
+     1+hh**2/4.d0*(5.d0*sd(i-1)+2.d0*sd(i-2)+5.d0*sd(i-3))
+  10  ukl(i)=2.d0*ukl(i-1)-ukl(i-2)+
+     1hh**2/12.d0*(sd(i)+10.d0*sd(i-1)+sd(i-2))
+c set initial values for the next block
+      do 20 j=2,jblock
+      hh=2.d0*hh
+      jj=(j-1)*nbl
+      y(jj+1)=ukl(jj+1)
+      y(jj)=ukl(jj-1)
+      y(jj-1)=ukl(jj-3)
+      y(jj-2)=ukl(jj-5)
+      do 30 jk=2,nb1
+      i=jj+jk
+      y(i)=y(i-1)+y(i-3)-y(i-4)+
+     1hh**2/4.d0*(5.d0*sdd(i-1)+2.d0*sdd(i-2)+5.d0*sdd(i-3))
+      y(i)=2.d0*y(i-1)-y(i-2)+
+     1hh**2/12.d0*(sdd(i)+10.d0*sdd(i-1)+sdd(i-2))
+  30  ukl(i)=y(i)
+  20  continue
+      return
+      end