subroutine starck(rcap,rsmall,thetad,vh3)
        implicit double precision(a-h,o-z)
        dimension p(5)
        data ap1/0.18734d0/,ap2/3.79103d0/,ap3/0.67432d0/
        data ap4/2.41997d0/
        data ap5/0.47723d0/,ap6/0.09074d0/,ap7/0.05819d0/
        data ap8/-2.35080d0/
        data ap9/-1.33464d0/,ap10/-0.06787d0/,ap11/-0.55760d0/
        data ap12/1.97391d0/,ap13/0.45551d0/
        data c1/35.42134d-3/,c2/-5.16557d-3/,c3/8.94017d-3/
        data c4/9.58594d-3/
        data c5/-4.42210d-3/,c6/3.11961d-3/,c7/-1.88474d-3/
        data c8/6.34267d-3/
        data c9/15.40166d-3/,c10/4.90275d-3/

        pi=dacos(-1.d0)
c
c       long range potential
c       vlr= sum over nl cnl(smallr)
c       c32, c40, c42, c54, c60, c62 as a function of small r
        theta=pi*thetad/180.d0
        r1f=rsmall
        r2f=dsqrt((r1f/2.d0)**2+rcap**2-2.0d0*rcap*(r1f/2.d0)
     1  *dcos(theta))
        r3f=dsqrt((r1f/2.d0)**2 +rcap**2 -2.0d0*rcap*(r1f/2.d0)*
     1  dcos(pi-theta))
c
        diff=r2f-r1f
        epart=ap1*dexp(-ap2*(diff/2.d0)**2)
        t1=r2f+r1f
        t2=dsqrt(diff**2+4.d0*epart**2)
        r1s=(t1-t2)/2.d0
        r2s=(t1+t2)/2.d0
        r3s=r3f
c
        sr=r1s
        cr=dsqrt((r2s**2+r3s**2-2.d0*(sr/2.d0)**2)/2.d0)
c       y=((sr/2.d0)**2+cr**2-r2s**2)/(2.d0*(sr/2.d0)*cr)
        an=(sr/2.d0)**2+cr**2-r2s**2
        ad=2.d0*(sr/2.d0)*cr
        y=an/ad
c       write(*,*)sr,cr,gama,sr,cr,r2s
c       do 10 i=1,20
c
        srlr1=sr-1.4d0
c       srlr1=rsmall-1.4d0
        srlr2=srlr1*srlr1
        srlr3=srlr2*srlr1
        srlr4=srlr3*srlr1
        srlr5=srlr4*srlr1
c       evaluation of c32(sr) etc
        c32sr=0.4555d0 + (0.5299d0*srlr1) + (0.0443d0*srlr2
     1  )+(-0.0909d0*srlr3)
     2  +(-0.0203d0*srlr4) + (0.0090d0*srlr5)
        c40sr=2.5897d0+ (2.1743d0*srlr1)+(0.4554d0*srlr2)
     1  +(-0.1678d0*srlr3)
     2  + (-0.1662d0*srlr4) + (0.0425d0*srlr5)
        c42sr=0.6053d0 + (1.1484d0*srlr1)+(0.6930d0*srlr2)
     1  +(-0.0573d0*srlr3)
     2  +(-0.2289d0*srlr4)+(0.0470d0*srlr5)
        c54sr=0.2745d0+(0.7337d0*srlr1)+(0.7138d0*srlr2)
     1  +(0.1425d0*srlr3)
     2  +(-0.1774d0*srlr4)+(0.0176d0*srlr5)
        c60sr=90.6083d0+(67.5401d0*srlr1)+(2.7573d0*srlr2)
     1  +(25.3497d0*srlr3)
     2  +(-24.3979d0*srlr4)+(5.0703d0*srlr5)
        c62sr=23.1469d0+(44.7755d0*srlr1)+(30.4350d0*srlr2)
     1  +(6.7026d0*srlr3)
     2  +(-1.6864d0*srlr4)+(-0.4605d0*srlr5)
c evaluation of dn"s. first convert capr and smallr into cr1 and sr1
c       read(5,*)ap1,ap2,ap3,ap4,ap5,ap6,ap7,ap8,ap9,ap10,ap11,ap12,
c    1  ap13
c       read(5,*)c1,c2,c3,c4,c5,c6,c7,c8,c9,c10
c convert sr into s1
c       sreq=1.4d0
        sr1=1.d0-dexp(ap3*(1.4d0-sr))
        sr1=sr1/ap3
c convert cr into cr1
        cr1=cr-3.0d0
c evaluation of dn"s. we have n=3,4,5,and 6.
c evaluation of xn"s that is, dn(xn)
c evaluation of pl(cos(gama)). l=0,2,and 4.
c       y=dcos(gama)
        call pleg(5,y,p)
c evaluation of dn(xn)
c evaluation d3x3
        xt=cr*(ap4+(ap5*sr1)+(ap6*cr1)+(ap7*p(3)))
        x3=xt
        if(x3.lt.0.d0)then
        d3x3=0.0d0
        else
        d3x3=1.d0-dexp(-x3)*(1.d0+x3+(x3*x3/2.d0)+(x3*x3*x3/6.d0))
        endif
c evaluation d4x4
        x4=xt+(cr*ap8)
        if(x4.lt.0.d0)then
        d4x4=0.d0
        else
        x42=x4*x4
        x43=x42*x4
        x44=x43*x4
        d4x4=1.d0-dexp(-x4)*(1.d0+x4+(x42/2.d0)+(x43/6.d0)
     1  +(x44/24.d0))
        endif
c  evaluatiom of d5x5
        x5=xt+(cr*ap8*2.d0)
         if(x5.lt.0.d0)then
        d5x5=0.d0
        else
        x52=x5*x5
        x53=x52*x5
        x54=x53*x5
        x55=x54*x5
        d5x5=1.d0 -dexp(-x5)*(1.d0+x5+(x52/2.d0)+(x53/6.d0)
     1  +(x54/24.d0)+(x55/120.d0))
        endif
c evaluation of d6x6
        x6=xt+(cr*ap8*3.d0)
        if(x6.lt.0.d0)then
        d6x6=0.d0
        else
        x62=x6*x6
        x63=x62*x6
        x64=x63*x6
        x65=x64*x6
        x66=x63*x63
        d6x6=1.d0 - dexp(-x6)*(1.d0+x6+(x62/2.d0)+
     1  (x63/6.d0)+(x64/24.d0)+(x65/120.)+(x66/720.d0))
        endif
c evaluation of long range pot : vlr
        cr2=cr*cr
        cr3=cr2*cr
        cr4=cr3*cr
        cr5=cr4*cr
        cr6=cr3*cr3
        vc32=c32sr*d3x3*p(3)/cr3
        vc40=c40sr*d4x4*p(1)/cr4
        vc42=c42sr*d4x4*p(3)/cr4
        vc54=c54sr*d5x5*p(5)/cr5
        vc60=c60sr*d6x6*p(1)/cr6
        vc62=c62sr*d6x6*p(3)/cr6
        vlr=vc32+vc40+vc42+vc54+vc60+vc62
        vlr= -vlr
c evaluation for exchange repulsion
        aex1=c1+(c2*sr1)+(c3*p(3))+(c4*cr1)+(c5*sr1*sr1)
        aex=aex1+(c6*p(3)*cr1)+(c7*cr1*cr1)+(c8*sr1*p(3))
        exarg=cr1*(ap9+ap10*p(3))
        vex=aex*dexp(exarg)
c evaluation for chem binding (resonance) energy
        rsarg=ap11*cr1+ap12*sr1
        betars=dexp(rsarg)*ap13*y
       term1=1.d0-dsqrt(1.d0+4.d0*betars*betars)
        term2=c9+c10*y*y
        vrs=term1*term2
c evaluation for asymptotic energy of h2
        san=1.-dexp(-1.5d0*((sr/1.4d0)-1.d0))
        san=san/1.5d0
        vh2=-1.17364d0 + (-0.00163d0*san) + (0.36017d0*san**2)
     1  +(-0.02754d0*san**3)+(0.15075d0*san**4)+(0.01849d0*san**5)
        vh2=vh2+1.17364d0
c evaluation of total interaction energy
         vas=vh2
         vt=vlr+vex+vrs+vas
         vh3=vt
         return
         end

         subroutine pleg(lmax,y,p)
         implicit double precision(a-h,o-z)
         dimension p(5)
         p(1)=1.d0
         p(2)=y
         if(lmax.le.1)return
         do 10 l=2,lmax
         lp=l+1
         lm=l-1
         xl=l
         p(lp)=((lm+l)*y*p(l)-lm*p(lm))/xl
 10      continue
         return
         end