c This file containes the following subroutines:
c ----------------------------------------------
c ldebas(ri,m,n,dbl,vl) : ri(m),dbl(6),vl(6): gives 6 type-n DbasL and Vlond.
c debas(ri,m,n,dbs) : ri(m),dbs(m): gives the m dEbas(H2) for vh2k and vh2s.
c getdeb(r,e,n,ypn,ep,s,sp) : r(n),e(n),ep(n),s(n),sp(n): initialize debas.
c spline(x,y,n,yp1,ypn,y2) : x(n),y(n),y2(n): sets up the spline array y2.
c splint(xa,ya,y2a,n,x,y) : xa(n),ya(n),y2a(n): gets value y of spline at x.
c splind(xa,ya,y2a,n,x,y,yp,ypp) : same as splint, but gets derivatives too.
c geth2n : initializes spline for vh2k, the Kolos et al H2 potential curve.
c vh2e3n(r,vh2k,vh2s,e30,e31,e32) : returns 2 singlet and 3 triplet values.
c vh2e3d(r,vh2k,vh2s,e30,e31,e32,vkp,vsp,vkpp,vspp) : same, with derivatives.
c---+*---|----+----|----+----|----+----|----+----|----+----|----+----|----+----
c
      subroutine ldebas(ri,m,n,dbl,vl)
c     --------------------------------
c
c Note arrays are of size: ri(m), dbl(6), vl(6);  m = 3 (for H3) or 6 (for H4).
c
c dBl = dVLondon(basis) ,  vl = VLondon(true)   for H3 or H4,
c
c where the triplet basis correction is here set to zero (ignored).
c
c Note that for H4 the six distances must be in the order A,B,C,D,E,F, i.e.,
c in the order r12,r23,r34,r24,r14,r13 for distances between atoms 1,2,3,4.
c
c dbi = dEbas(H2) ,   z = dEbas(3E) [triplet]
c
c Note 6 values are returned for "dbl" and "vl": three using vh2k with the
c respective triplets e30, e31, e32; three using vh2s and these triplets.
c vh2k is the (new) spline H2 curve, while vh2s is the Schwenke H2 curve;
c e30 is the Truhlar & Horowitz triplet, e31 is the corrected lin-sym triplet
c (using the 25 linsym points with Schwenke basis correction) used for the
c old DbasL, and e32 is the final linsym+3gengenbach triplet used for the new
c (final) DbasL.
c
c Note that for H3 and H4 basis correction, vh2s and e32 were used, i.e.,
c the values returned in dbl(6) and vl(6).
c
c n=-3 for Partridge's 6s/6p/4d/2f basis, n=-2 for Partridge's 6s/5p/3d/1f,
c n=-1 for Liu's 4s/3p/2d? Slater basis, n=0 for our (Gaussian) 4s/3p/2d basis,
c n=1 for 4s/3p/1d, n=2 for 4s/3p, n=3 for Schwenke's 4s/2p/1d basis
c
      implicit real*8(a-h,o-z)
      dimension ri(m),dbi(12),dbs(6),s(6,6),t(6,6),dbl(6),vl(6),
     $     e3(0:2),vh2(2)
      equivalence (vh2(1),vh2k), (vh2(2),vh2s), (dbi(7),dbs(1))
c
      if(m.ne.3.and.m.ne.6)stop 'STOP -- LDEBAS: Neither H3 nor H4'
      do i=1,m
         dbi(i)=ri(i)
      enddo
      call debas(dbi,m,n,dbs)
      do k=1,6
         dbl(k)=0.d0
         vl(k)=0.d0
      enddo
      do i=1,m
         call vh2e3n(ri(i),vh2(1),vh2(2),e3(0),e3(1),e3(2))
         do k=1,6
            jv=(k+2)/3
            je=k+2-jv*3
            ji=i-6+jv*6
            s(i,k)=(vh2(jv)-e3(je))*.5d0
            t(i,k)=(vh2(jv)-dbi(ji)-e3(je))*.5d0
            dbl(k)=dbl(k)+dbi(ji)
            vl(k)=vl(k)+vh2(jv)+e3(je)
         enddo
      enddo
      do k=1,6
         if(m.eq.3)then
            tau=s(2,k)-s(1,k)
            pi=s(3,k)-s(1,k)
            sig=s(3,k)-s(2,k)
            tt=t(2,k)-t(1,k)
            pp=t(3,k)-t(1,k)
            ss=t(3,k)-t(2,k)
         else
            tau=s(1,k)+s(3,k)-s(6,k)-s(4,k)
            pi=s(6,k)+s(4,k)-s(5,k)-s(2,k)
            sig=s(5,k)+s(2,k)-s(1,k)-s(3,k)
            tt=t(1,k)+t(3,k)-t(6,k)-t(4,k)
            pp=t(6,k)+t(4,k)-t(5,k)-t(2,k)
            ss=t(5,k)+t(2,k)-t(1,k)-t(3,k)
         endif
         sj=dsqrt((tau**2+pi**2+sig**2)*.5d0)
         dbl(k)= ( dsqrt((tt**2+pp**2+ss**2)*.5d0) - sj ) + dbl(k)*.5d0
         vl(k)= vl(k)*.5d0 - sj
      enddo
      return
      end
c
c---+----|----+----|----+----|----+----|----+----|----+----|----+----|----+----
c
      subroutine debas(ri,m,n,dbs)
c     ----------------------------
c
c RI(M) is an array of M pair distances; N=-3 for Partridge's 6s6p4d2f basis,
c N=-2 for Partridge's 6s5p3d1f basis, N=-1 for Liu's 4s/3p/2d Slater basis,
c N=0 for our (Gaussian) 4s/3p/2d basis, N=1 for 4s/3p/1d, N=2 for 4s/3p,
c N=3 for Schwenke's 4s/2p/1d basis;
c dEbas values relative to the Kolos et al H2 potential curve are returned
c in RI, replacing their distance values; the dEbas values relative to the
c Schwenke H2 potential curve are returned in DBS, another array of size M.
c  Note that the arrays RI(M),DBS(M) must be REAL*8 (i.e. double precision).
c
      implicit real*8 (a-h,o-z)
c
      parameter (nptsp=40, nptspd=40, nptd2=40, nptsla=18, nptsch=21,
     $     schmul=16.d0/15.d0,
c-dontuse     $     np6531=34,
     $     np6531=33, np6642=16, nlo=-3, nhi=3, numn=nhi-nlo+1)
      common/isinit/init(nlo:nhi)
      common/dbmat/rvsp(nptsp),rvspd(nptspd),despd(nptspd),desp(nptsp),
     $     rvd2(nptd2),ded2(nptd2),depspd(nptspd),depsp(nptsp),
     $     depd2(nptd2),rvslat(nptsla),deslat(nptsla),depsla(nptsla),
     $     rvsch(nptsch),desch(nptsch),depsch(nptsch),seslat(nptsla),
     $     sepsla(nptsla),sed2(nptd2),sepd2(nptd2),sespd(nptspd),
     $     sepspd(nptspd),sesp(nptsp),sepsp(nptsp),sesch(nptsch),
     $     sepsch(nptsch),rv6642(np6642),rv6531(np6531),de6642(np6642),
     $     de6531(np6531),dep6642(np6642),dep6531(np6531),
     $     se6642(np6642),se6531(np6531),sep6642(np6642),
     $     sep6531(np6531),espd(nptspd),esp(nptsp),ed2(nptd2),
     $     esch(nptsch),eslat(nptsla),e6531(np6531),e6642(np6642)
      dimension ri(m),dbs(m)
c
      data init/numn*0/
c-old;      data rvspd/.5d0,.6d0,.8d0,1.d0,1.2d0,1.35d0,1.4d0,1.45d0,1.6d0,
c-old;     $     1.8d0,2.d0,2.25d0,2.5d0,2.75d0,3.d0,3.25d0,3.5d0,3.75d0,
c-old;     $     4.d0,4.5d0,5.d0,6.d0,7.d0/
c-old;      data despd/3.19482d0,2.57879d0,1.82954d0,1.47481d0,1.27697d0,
c-old;     $     1.18678d0,1.16652d0,1.15108d0,1.12962d0,1.11156d0,1.08749d0,
c-old;     $     1.08960d0,1.11917d0,1.15180d0,1.16298d0,1.13277d0,1.05281d0,
c-old;     $     .93372d0,.79619d0,.53429d0,.33976d0,.13536d0,.05552d0/
      data rvspd/.5d0,.55d0,.6d0,.65d0,.7d0,.75d0,.8d0,.9d0,1.d0,
     $     1.1d0,1.2d0,1.3d0,1.35d0,1.4d0,1.45d0,1.5d0,1.6d0,1.7d0,
     $     1.8d0,1.9d0,2.d0,2.2d0,2.4d0,2.6d0,2.8d0,3.d0,3.2d0,3.4d0,
     $     3.6d0,3.8d0,4.d0,4.2d0,4.4d0,4.6d0,4.8d0,5.d0,5.5d0,6.d0,
     $     6.5d0,7.d0/
      data espd/.47655654d0,.34009549d0,.23294378d0,.14797388d0,
     $     .08010659d0,.02562994d0,-.01822677d0,-.08202029d0,
     $     -.12306458d0,-.14869388d0,-.16365799d0,-.17113499d0,
     $     -.17277670d0,-.17330896d0,-.17290577d0,-.17171460d0,
     $     -.16745356d0,-.16133526d0,-.15395702d0,-.14575297d0,
     $     -.13704531d0,-.11904572d0,-.10131698d0,-.08465779d0,
     $     -.06952662d0,-.05616311d0,-.04465662d0,-.03498464d0,
     $     -.02703726d0,-.02064249d0,-.01559363d0,-.01167375d0,
     $     -.00867383d0,-.00640545d0,-.00470695d0,-.00344543d0,
     $     -.00156047d0,-.00070012d0,-.00031407d0,-.00014232d0/
c
c-old;      data rvsp/.5d0,.6d0,.8d0,1.d0,1.2d0,1.35d0,1.4d0,1.45d0,1.6d0,
c-old;     $     2.d0,2.5d0,3.d0,3.5d0,4.d0,5.d0,6.d0,7.d0/
c-old;      data desp/ 6.51816d0,4.97367d0,3.41528d0,2.84223d0,2.66682d0,
c-old;     $     2.64318d0,2.64528d0,2.65010d0,2.66903d0,2.62630d0,2.37778d0,
c-old;     $     1.99641d0,1.49847d0,1.00110d0,.37639d0,.14194d0,.05677d0/
      data rvsp/.5d0,.55d0,.6d0,.65d0,.7d0,.75d0,.8d0,.9d0,1.d0,
     $     1.1d0,1.2d0,1.3d0,1.35d0,1.4d0,1.45d0,1.5d0,1.6d0,1.7d0,
     $     1.8d0,1.9d0,2.d0,2.2d0,2.4d0,2.6d0,2.8d0,3.d0,3.2d0,3.4d0,
     $     3.6d0,3.8d0,4.d0,4.2d0,4.4d0,4.6d0,4.8d0,5.d0,5.5d0,6.d0,
     $     6.5d0,7.d0/
      data esp/.47988053d0,.34288828d0,.23533910d0,.15007267d0,
     $     .08198338d0,.02734011d0,-.01664088d0,-.08058873d0,
     $     -.12169715d0,-.14733208d0,-.16226818d0,-.16970183d0,
     $     -.17132034d0,-.17183024d0,-.17140679d0,-.17019848d0,
     $     -.16591417d0,-.15978324d0,-.15239791d0,-.14419658d0,
     $     -.13550649d0,-.11758659d0,-.09998356d0,-.08347939d0,
     $     -.06851848d0,-.05532965d0,-.04399187d0,-.03447180d0,
     $     -.02665230d0,-.02035941d0,-.01538872d0,-.01152697d0,
     $     -.00856954d0,-.00633171d0,-.00465497d0,-.00340879d0,
     $     -.00154502d0,-.00069359d0,-.00031117d0,-.00014100d0/
c
c-old;      data rvd2/.5d0,.6d0,.8d0,1.d0,1.2d0,1.35d0,1.4d0,1.45d0,1.6d0,
c-old;     $     1.8d0,2.d0,2.2d0,2.5d0,3.d0,3.5d0,4.d0,4.5d0,5.d0,6.d0,7.d0/
c-old;      data ded2/1.94439d0,1.63641d0,1.15683d0,.90711d0,.78679d0,
c-old;     $     .73588d0,.72359d0,.71326d0,.69272d0,.68460d0,.68931d0,
c-old;     $     .69757d0,.70459d0,.69146d0,.63550d0,.52700d0,.39154d0,
c-old;     $     .27053d0,.12104d0,.05257d0/
      data rvd2/.5d0,.55d0,.6d0,.65d0,.7d0,.75d0,.8d0,.9d0,1.d0,
     $     1.1d0,1.2d0,1.3d0,1.35d0,1.4d0,1.45d0,1.5d0,1.6d0,1.7d0,
     $     1.8d0,1.9d0,2.d0,2.2d0,2.4d0,2.6d0,2.8d0,3.d0,3.2d0,3.4d0,
     $     3.6d0,3.8d0,4.d0,4.2d0,4.4d0,4.6d0,4.8d0,5.d0,5.5d0,6.d0,
     $     6.5d0,7.d0/
      data ed2/.47530600d0,.33902040d0,.23200134d0,.14713106d0,
     $     .07933778d0,.02491594d0,-.01889953d0,-.08263397d0,
     $     -.12363237d0,-.14922043d0,-.16414826d0,-.17159665d0,
     $     -.17322768d0,-.17375196d0,-.17334365d0,-.17215016d0,
     $     -.16789050d0,-.16177217d0,-.15438400d0,-.14616458d0,
     $     -.13744351d0,-.11943442d0,-.10171911d0,-.08508656d0,
     $     -.06998228d0,-.05663465d0,-.04512416d0,-.03542419d0,
     $     -.02742827d0,-.02097352d0,-.01586286d0,-.01188604d0,
     $     -.00883746d0,-.00652946d0,-.00480011d0,-.00351463d0,
     $     -.00159250d0,-.00071438d0,-.00032035d0,-.00014528d0/
c
      data rvslat/.6d0,.8d0,1.d0,1.2d0,1.4d0,1.6d0,1.8d0,2.d0,2.2d0,
     $     2.4d0,2.6d0,2.8d0,3.d0,3.2d0,3.5d0,4.d0,5.d0,6.d0/
      data deslat/.975d0,.897d0,.841d0,.795d0,.763d0,.741d0,.740d0,
     $     .738d0,.726d0,.715d0,.697d0,.673d0,.644d0,.602d0,.540d0,
     $     .425d0,.216d0,.072d0/
c
      data rvsch/.5d0,.6d0,.8d0,1.d0,1.2d0,1.4d0,1.6d0,1.8d0,2.d0,
     $     2.25d0,2.5d0,2.75d0,3.d0,3.25d0,3.5d0,3.75d0,4.d0,4.5d0,
     $     5.d0,6.d0,7.d0/
      data desch/5.10597d0,4.05947d0,2.68009d0,1.93915d0,1.56775d0,
     $     1.39229d0,1.32078d0,1.31002d0,1.33892d0,1.40877d0,
     $     1.48961d0,1.55167d0,1.56683d0,1.51837d0,1.40833d0,
     $     1.25454d0,1.08056d0,.74730d0,.48916d0,.19507d0,.07950d0/
c
c-old;      data rv6531/1.d0,1.2d0,1.4d0,1.449d0,1.6d0,1.7d0,1.8d0,2.d0,
c-old;     $     2.2d0,2.6d0,3.d0,4.d0,6.d0,8.d0/
c-old;      data de6531/.65914d0,.55307d0,.48113d0,.46742d0,.43244d0,
c-old;     $     .41418d0,.39897d0,.37515d0,.35675d0,.32684d0,.29780d0,
c-old;     $     .17705d0,.03030d0,.00607d0/
c-dontuse;      data rv6531/.8d0,
      data rv6531/.9d0,.95d0,1.d0,1.05d0,1.1d0,1.15d0,1.2d0,
     $     1.25d0,1.3d0,1.4d0,1.449d0,1.5d0,1.6d0,1.7d0,1.75d0,1.8d0,
     $     1.9d0,2.d0,2.2d0,2.4d0,2.6d0,2.8d0,3.d0,3.4d0,3.8d0,4.d0,
     $     4.5d0,5.d0,5.5d0,6.d0,6.5d0,7.d0,8.d0/
c-dontuse;      data e6531/-.01964331d0,  ! gives dE=.41329 mh, should be .77
      data e6531/-.08290665d0,-.10569977d0,-.12388052d0,
     $     -.13825085d0,-.14945650d0,-.15802352d0,-.16438212d0,
     $     -.16888723d0,-.17183355d0,-.17399454d0,-.17360611d0,
     $     -.17240057d0,-.16815089d0,-.16204450d0,-.15848532d0,
     $     -.15466972d0,-.14646449d0,-.13775776d0,-.11977532d0,
     $     -.10208144d0,-.08546434d0,-.07037029d0,-.05702837d0,
     $     -.03581724d0,-.02134469d0,-.01621290d0,-.00787702d0,
     $     -.00371123d0,-.00172715d0,-.00080524d0,-.00038078d0,
     $     -.00018488d0,-.00004951d0/
c
c-old;      data rv6642/.8d0,1.d0,1.2d0,1.4d0,1.449d0,1.6d0,1.8d0,2.d0,
c-old;     $     2.4d0,3.d0,4.d0,6.d0,8.d0/
c-olderr;      data de6642/.47162d0,.41953d0,.37905d0,.34248d0,.33527d0,
c-olderr;     $     .31423d0,.28849d0,.26712d0,.23317d0,.18480d0,.11061d0,
c-olderr;     $     .02100d0,.00438d0/
      data rv6642/.8d0,1.d0,1.2d0,1.4d0,1.449d0,1.6d0,1.75d0,1.8d0,
     $     2.d0,2.2d0,2.4d0,2.6d0,3.d0,4.d0,6.d0,8.d0/
      data e6642/-.01964589d0,-.12412013d0,-.16455614d0,-.17416126d0,
     $     -.17376375d0,-.16828665d0,-.15860666d0,-.15478735d0,
     $     -.13786438d0,-.11987628d0,-.10218103d0,-.08556502d0,
     $     -.05714138d0,-.01627934d0,-.00081273d0,-.00005118d0/
c
      if(n.lt.nlo.or.n.gt.nhi)stop ' STOP -- DEBAS n-value error'
      dhi=7.d0
      if(n.lt.0)then
         dhi=6.d0
         if(n.lt.-1)dhi=8.d0
      endif
      do i=1,m
         rv=ri(i)
         if(rv.lt.0.3d0)stop ' STOP -- Error: value of R < 0.3'
         if(rv.gt.dhi)then
            call vh2e3n(rv,vh2,vh2s,e30,e31,e32)
            if(n.lt.-1)then
               rvm=rv
               if(n.eq.-3)rvm=rv-.6d0
               d=vh2/(dexp((11.5d0-rvm)*.599673d0)+1.d0)
            else if(n.eq.-1)then
               d=vh2*dmin1(.026d0*rv-.0698d0,1.d0)
            else if(n.eq.3)then
c-schmul=1               d=vh2/(dexp((7.51d0-rv)*.780175d0)+1.d0)
               d=vh2/(dexp((7.355d0-rv)*.81d0)+1.d0)
            else
               if(n.eq.0)then
                  rvm=rv-dexp((3.16d0-rv)*.57498d0)
               else if(n.eq.2)then
                  rvm=rv+dexp((2.54d0-rv)*.69325d0)
               else
                  rvm=rv
               endif
               d=vh2/(dexp((8.377d0-rvm)*.68364d0)+1.d0)
            endif
            s=d+(vh2s-vh2)
         else if(n.lt.-2)then
            if(init(-3).eq.0)then
               ypn=-.002635d0
               call getdeb(rv6642,de6642,np6642,ypn,dep6642,se6642,
     $              sep6642,e6642,n)
               init(-3)=1
            endif
            call splint(rv6642,de6642,dep6642,np6642,rv,d)
            d=-1.d-3*d
            call splint(rv6642,se6642,sep6642,np6642,rv,s)
            s=-1.d-3*s
         else if(n.eq.-2)then
            if(init(-2).eq.0)then
               ypn=-.003762d0
               call getdeb(rv6531,de6531,np6531,ypn,dep6531,se6531,
     $              sep6531,e6531,n)
               init(-2)=1
            endif
            call splint(rv6531,de6531,dep6531,np6531,rv,d)
            d=-1.d-3*d
            call splint(rv6531,se6531,sep6531,np6531,rv,s)
            s=-1.d-3*s
         else if(n.eq.-1)then
            if(init(-1).eq.0)then
               ypn=-.0858d0
               call getdeb(rvslat,deslat,nptsla,ypn,depsla,seslat,
     $              sepsla,eslat,n)
               init(-1)=1
            endif
            call splint(rvslat,deslat,depsla,nptsla,rv,d)
            d=-1.d-3*d
            call splint(rvslat,seslat,sepsla,nptsla,rv,s)
            s=-1.d-3*s
         else if(n.eq.0)then
            if(init(0).eq.0)then
               ypn=-.043972d0
               call getdeb(rvd2,ded2,nptd2,ypn,depd2,sed2,sepd2,ed2,n)
               init(0)=1
            endif
            call splint(rvd2,ded2,depd2,nptd2,rv,d)
            d=-1.d-3*d
            call splint(rvd2,sed2,sepd2,nptd2,rv,s)
            s=-1.d-3*s
         else if(n.eq.1)then
            if(init(1).eq.0)then
               ypn=-.048769d0
               call getdeb(rvspd,despd,nptspd,ypn,depspd,sespd,sepspd,
     $              espd,n)
               init(1)=1
            endif
            call splint(rvspd,despd,depspd,nptspd,rv,d)
            d=-1.d-3*d
            call splint(rvspd,sespd,sepspd,nptspd,rv,s)
            s=-1.d-3*s
         else if(n.eq.2)then
            if(init(2).eq.0)then
               ypn=-.050982d0
               call getdeb(rvsp,desp,nptsp,ypn,depsp,sesp,sepsp,esp,n)
               init(2)=1
            endif
            call splint(rvsp,desp,depsp,nptsp,rv,d)
            d=-1.d-3*d
            call splint(rvsp,sesp,sepsp,nptsp,rv,s)
            s=-1.d-3*s
         else if(n.eq.3)then
            if(init(3).eq.0)then
               do j=1,nptsch
                  desch(j)=desch(j)*schmul
               enddo
c-schmul=1               ypn=-.071828d0
               ypn=-.0770175d0
               call getdeb(rvsch,desch,nptsch,ypn,depsch,sesch,sepsch,
     $              esch,n)
               init(3)=1
            endif
            call splint(rvsch,desch,depsch,nptsch,rv,d)
            d=-1.d-3*d
            call splint(rvsch,sesch,sepsch,nptsch,rv,s)
            s=-1.d-3*s
         endif
         ri(i)=d
         dbs(i)=s
      enddo
      return
      end
c
      subroutine getdeb(r,d,np,ypn,ep,s,sp,e,n)
c     -----------------------------------------
c
      implicit real*8(a-h,o-z)
      dimension r(np),d(np),ep(np),s(np),sp(np),e(np)
c
      if(n.eq.-1.or.n.eq.3)then
         call spline(r,d,np,1.d37,ypn,ep)
         do i=1,np
            call vh2e3d(r(i),vh2k,vh2s,e30,e31,e32,vkp,vsp,vkpp,vspp)
            e(i)=vh2k+1.d-3*d(i)
            s(i)=(vh2k-vh2s)*1.d3+d(i)
         enddo
      else
         do i=1,np
            call vh2e3d(r(i),vh2k,vh2s,e30,e31,e32,vkp,vsp,vkpp,vspp)
            d(i)=(e(i)-vh2k)*1.d3
            s(i)=(e(i)-vh2s)*1.d3
         enddo
         call spline(r,d,np,1.d37,ypn,ep)
      endif
      yps=(vkp-vsp)*1.d3+ypn
      call spline(r,s,np,1.d37,yps,sp)
      return
      end
c
c---+*---|----+----|----+----|----+----|----+----|----+----|----+----|----+----
c
      subroutine splint(xa,ya,y2a,n,x,y)
c     ----------------------------------
c
      implicit real*8(a-h,o-z)
      dimension xa(n),ya(n),y2a(n)
c
      call splind(xa,ya,y2a,n,x,y,yp,ypp)
      return
      end
c
      subroutine splind(xa,ya,y2a,n,x,y,yp,ypp)
c     -----------------------------------------
c
c Double precision spline evaluation program: given arrays xa(n), ya(n),
c and second derivative array y2a(n), returns the cubic-spline interpolated
c value y at the given value x.
c
      implicit real*8(a-h,o-z)
      dimension xa(n),ya(n),y2a(n)
      klo=1
      khi=n
      do while(khi-klo.gt.1)
         k=(khi+klo)/2
         if(xa(k).gt.x)then
	    khi=k
         else
	    klo=k
         endif
      enddo
      h=xa(khi)-xa(klo)
      if(h.le.0.d0)stop 'STOP -- SPLINT Error: x-values not in order'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      cah=(a**2-1.d0)*h/6.d0
      dbh=(b**2-1.d0)*h/6.d0
      y= a*ya(klo) + b*ya(khi)
     $     + (a*cah*y2a(klo)+b*dbh*y2a(khi))*h
      yp= (ya(khi)-ya(klo))/h - (3.d0*cah+h/3.d0)*y2a(klo)
     $     + (3.d0*dbh+h/3.d0)*y2a(khi)
      ypp=a*y2a(klo)+b*y2a(khi)
      return
      end
c
      subroutine spline(x,y,n,yp1,ypn,y2)
c     -----------------------------------
c
c Double precision spline creation program: give arrays x(n) and y(n),
c returns array y2(n) of second derivatives for cubic spline.  Must be
c given values yp1 and ypn of the first derivative at the endpoints; if
c greater than 1.E35, then "natural spline" is used, i.e. second derivative
c zero at the endpoint(s).  Size n of array must not exceed nmax=100.
c
      implicit real*8(a-h,o-z)
      parameter (nmax=100)
      dimension x(n),y(n),y2(n),u(nmax)
      if(n.gt.nmax)stop 'STOP -- Error: array size too big'
      if(yp1.gt.0.999999d35)then
         y2(1)=0.d0
         u(1)=0.d0
      else
         y2(1)=-0.5d0
         u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do i=2,n-1
         sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
         if(sig.le.0.d0)
     $        stop 'STOP -- Error: X-values not in order'
         p=sig*y2(i-1)+2.d0
         y2(i)=(sig-1.d0)/p
         u(i)= ( 6.d0 * ( (y(i+1)-y(i))/(x(i+1)-x(i))
     $        - (y(i)-y(i-1))/(x(i)-x(i-1)) ) / (x(i+1)-x(i-1))
     $        - sig*u(i-1) ) / p
      enddo
      if(x(n).le.x(n-1))stop 'STOP -- Error: X-values not in order'
      if(ypn.gt.0.999999d35)then
         qn=0.d0
         un=0.d0
      else
         qn=0.5d0
         un=3.d0/(x(n)-x(n-1))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
      do i=n-1,1,-1
         y2(i)=y2(i)*y2(i+1)+u(i)
      enddo
      return
      end
c
c---+*---|----+----|----+----|----+----|----+----|----+----|----+----|----+----
c
      subroutine vh2e3n(r,vh2k,vh2s,e30,e31,e32)
c     ------------------------------------------
c
      implicit real*8(a-h,o-z)
c
      call vh2e3d(r,vh2k,vh2s,e30,e31,e32,vkp,vsp,vkpp,vspp)
      return
      end
c
      subroutine vh2e3d(r,vh2k,vh2s,e30,e31,e32,vkp,vsp,vkpp,vspp)
c     ------------------------------------------------------------
c
c For separation r, this subroutine returns the Kolos et al H2 potential vh2k,
c the Schwenke H2 potential vh2s, the Truhlar & Horowitz triplet e30, the old
c corrected lin-sym triplet e31 (using the 25 linsym points with Schwenke
c basis correction), and the final linsym+3gengenbach triplet e32.
c
      implicit real*8(a-h,o-z)
      parameter ( r0=3.5284882d0, DD=0.160979391d0, c6=6.499027d0,
     $     c8=124.3991d0, c10=3285.828d0, r02=r0*r0, r04=r02*r02,
     $     r06=r04*r02 , v6=-9.163021333d0, v8=642.104267d0,
     $     v10=-60100.21337d0 )
      common/h2valn/rn(83),en(83),enp(83)
      common/inspln/initsp
      dimension a(17),rpow(0:15)
      data initsp/0/
      data a/0.03537359271649620d0,  2.013977588700072d0,
     $      -2.827452449964767d0,    2.713257715593500d0,
     $      -2.792039234205731d0,    2.166542078766724d0,
     $      -1.272679684173909d0,    0.5630423099212294d0,
     $      -0.1879397372273814d0,   0.04719891893374140d0,
     $      -0.008851622656489644d0, 0.001224998776243630d0,
     $      -1.227820520228028d-04,  8.638783190083473d-06,
     $      -4.036967926499151d-07,  1.123286608335365d-08,
     $      -1.406619156782167d-10/
c
      if(r.le.0.d0)stop 'STOP -- Error in R value'
      alphar=a(1)/r+a(2)+a(3)*r
      dalphr=a(3)
      ddalph=0.d0
      rpow(0)=1.d0
      rpow(1)=r
      do i=2,15
         rpow(i)=rpow(i-1)*r
         alphar=alphar+a(i+2)*rpow(i)
         dalphr=dalphr+i*a(i+2)*rpow(i-1)
         ddalph=ddalph+i*(i-1)*a(i+2)*rpow(i-2)
      enddo
      dalphr=dalphr-a(1)/rpow(2)
      ddalph=ddalph+2.d0*a(1)/rpow(3)
      exalph=dexp(alphar)
      rr6=rpow(6)+r06
      rr4=rpow(4)+r04
      rr2=rpow(2)+r02
      vsr = DD*( exalph - 1.d0 )**2 - DD
      vlr = -C6/rr6 - C8/rr4**2 - C10/rr2**5
      vh2s=vsr+vlr
      dvsr = 2.d0 * DD*(exalph-1.d0) * exalph * dalphr
      dvlr = 6.d0*c6*rpow(5)/rr6**2 + 8.d0*c8*rpow(3)/rr4**3
     $     + 10.d0*c10*r/rr2**6
      vsp=dvsr+dvlr
      ddvsr = 2.d0 * DD * exalph * ( (2.d0*exalph-1.d0)*dalphr**2
     $     + (exalph-1.d0)*ddalph )
      ddvlr = c6 * ( 30.d0*rpow(4) - 72.d0*rpow(10)/rr6 ) / rr6**2
     $     + c8 * ( 24.d0*rpow(2) - 96.d0*rpow(6)/rr4 ) / rr4**3
     $     + c10 * ( 10.d0 - 120.d0*rpow(2)/rr2 ) / rr2**6
      vspp=ddvsr+ddvlr
c
      if(r.le.10.d0)then
         if(initsp.eq.0)then
            initsp=1
            call geth2n
         endif
         call splind(rn,en,enp,83,r,vh2k,vkp,vkpp)
         vh2k=vh2k+1.d0
      else
         vh2k = v6/rpow(6) + v8/rpow(8) + v10/rpow(10)
         vkp = -6.d0*v6/rpow(7) - 8.d0*v8/rpow(9) - 10.d0*v10/rpow(11)
         vkpp = 42.d0*v6/rpow(8) + 72.d0*v8/rpow(10)
     $        + 121.d0*v10/rpow(12)
      endif
c
c Truhlar & Horowitz triplet:
c
      e30=-1.2148730613d0*((-1.46d0*r+1.d0)*r-1.5146634740d0)
     $     *dexp(-2.088442d0*r)
c
c all H3:
c
c-improved      e31=-1.27258594d0*( (-1.94132955d0+r*(1.d0-.6965810413d0*r))
c-improved     $     *dexp(-1.695978560d0*r) +.06028475976d0/r**1.899124608d0 )
c
c-nouse      e31=0.15980625d0*( (1.13293947d0+r*(1.d0-0.32753829d0*r))
c-nouse     $     *dexp(-0.74564871d0*r) +0.8992885717d0/r**2.15789966d0 )
c
c all H3 free:
c
c-nouse      e31=0.4304999809d0*( (0.1225373969d0+r*(1.d0-0.2176502848d0*r))
c-nouse     $     *dexp(-0.9826419772d0*r) +0.3022602970d0/r**2.932567027d0 )
c
c all H4 free:
c
c-nouse      e31=1.56632282d0*( (-1.11217629d0+r*(1.d0-0.28235740d0*r))
c-nouse     $     *dexp(-0.94521838d0*r) +0.33821022d0/r**1.60977620d0 )
c
c lin-sym triplet e31 (from 25 linsym pts: p.30 of Bill)
c
c      e31=-.281929d0*( (-3.063451d0+r*(1.d0-5.179881d0*r))
c     $     *dexp(-2.092891d0*r) -.16682303d0/r**3.0304533d0 )
c
c linsym: supposedly (but not actually) used for old DbasL:
c
      e31=-.50477840d0*( (-1.95782581d0+r*(1.d0-3.21921775d0*r))
     $     *dexp(-2.12139394d0*r) -.091947239d0/r**2.928044d0 )
c
c linsym: not used for old DbasL:
c
c      e31=-.50721542d0*( (-1.94818782d0+r*(1.d0-3.20929885d0*r))
c     $     *dexp(-2.12199549d0*r) -.091732080d0/r**2.925775d0 )
c
c linsym+3geng(?):
c
c-improved      e32=-0.59459228d0*( (-2.07517218d0+r*(1.d0-2.61134931d0*r))
c-improved     $     *dexp(-2.08300479d0*r) -.04351971938d0/r**3.51307678d0 )
c
c final linsym+3geng:
c
c-improved      e32=-0.66129429d0*( (-1.99434198d0+r*(1.d0-2.37604328d0*r))
c-improved     $     *dexp(-2.08107802d0*r) -0.0313032510d0/r**3.76546699d0 )
c
c final linsym+3geng:
c
c-improved      e32=-0.66053951d0*( (-1.99632184d0+r*(1.d0-2.37786221d0*r))
c-improved     $     *dexp(-2.08097071d0*r) -.0312982924309d0/r**3.76839356d0 )
c
c final linsym+3geng:   (rms=0.0474 unweighted and 0.0483 weighted).
c
      e32 = -0.0253496194d0*( -29.2302126444+r*(1.d0-50.7225015503d0*r)
     $     -12.2408908509d0/r**1.6733157383d0 )*dexp(-2.0452676876d0*r)
c
      return
      end
c
      subroutine geth2n
c     -----------------
c
c This initialization routine sets up the spline for computing vh2k.
c
      implicit real*8(a-h,o-z)
      common/h2valn/rn(83),en(83),enp(83)
c
      data rn/	0.4d0,	0.45d0,	0.5d0,	0.55d0,	0.6d0,	0.65d0,	0.7d0,
     $	0.75d0,	0.8d0,	0.9d0,	1.d0,	1.1d0,	1.2d0,	1.3d0,	1.35d0,
     $	1.4d0,	1.45d0,	1.5d0,	1.6d0,	1.7d0,	1.8d0,	1.9d0,	2.d0,
     $	2.1d0,	2.2d0,	2.3d0,	2.4d0,	2.5d0,	2.6d0,	2.7d0,	2.8d0,
     $	2.9d0,	3.d0,	3.1d0,	3.2d0,	3.3d0,	3.4d0,	3.5d0,	3.6d0,
     $	3.7d0,	3.8d0,	3.9d0,	4.d0,	4.1d0,	4.2d0,	4.3d0,	4.4d0,
     $	4.5d0,	4.6d0,	4.7d0,	4.8d0,	4.9d0,	5.d0,	5.1d0,	5.2d0,
     $	5.3d0,	5.4d0,	5.5d0,	5.6d0,	5.7d0,	5.8d0,	5.9d0,	6.d0,
     $	6.1d0,	6.2d0,	6.3d0,	6.4d0,	6.5d0,	6.6d0,	6.7d0,	6.8d0,
     $	6.9d0,	7.d0,	7.2d0,	7.4d0,	7.6d0,	7.8d0,	8.d0,	8.25d0,
     $	8.5d0,	9.d0,	9.5d0,	10.d0/
c
      data en/		-0.120230238d0,	-0.3509478d0,	-0.526638668d0,
     $	-0.6627745d0,	-0.769635349d0,	-0.8543625d0,	-0.9220271d0,
     $	-0.9763367d0,	-1.020056598d0,	-1.083643176d0,	-1.124539660d0,
     $	-1.150057312d0,	-1.164935191d0,	-1.172347100d0,	-1.173963679d0,
     $	-1.174475668d0,	-1.174057026d0,	-1.172855034d0,	-1.168583330d0,
     $	-1.162458684d0,	-1.155068695d0,	-1.1468506d0,	-1.138132913d0,
     $	-1.1291638d0,	-1.120132069d0,	-1.1111817d0,	-1.102422553d0,
     $	-1.0939381d0,	-1.085791176d0,	-1.0780284d0,	-1.070683159d0,
     $	-1.0637779d0,	-1.057326175d0,	-1.0513337d0,	-1.045799543d0,
     $	-1.0407172d0,	-1.036075238d0,	-1.0318579d0,	-1.028046110d0,
     $	-1.0246180d0,	-1.021549547d0,	-1.0188155d0,	-1.016389951d0,
     $	-1.0142466d0,	-1.012359608d0,	-1.0107041d0,	-1.009256107d0,
     $	-1.0079933d0,	-1.006894750d0,	-1.0059413d0,	-1.005115473d0,
     $	-1.0044015d0,	-1.003785212d0,	-1.0032539d0,	-1.002796446d0,
     $	-1.0024029d0,	-1.002064751d0,	-1.0017743d0,	-1.001524999d0,
     $	-1.0013111d0,	-1.001127673d0,	-1.0009704d0,	-1.000835537d0,
     $	-1.0007199d0,	-1.000620826d0,	-1.0005358d0,	-1.000462967d0,
     $	-1.0004005d0,	-1.000346789d0,	-1.0003007d0,	-1.000261141d0,
     $	-1.0002271d0,	-1.000197852d0,	-1.000150944d0,	-1.000116047d0,
     $	-1.000089970d0,	-1.000070382d0,	-1.000055582d0,	-1.0000420d0,
     $	-1.000032157d0,	-1.000019773d0,	-1.000012851d0,	-1.000008572d0/
c
      yp1=-5.305861654d0
      ypn=.000006371d0
      call spline(rn,en,83,yp1,ypn,enp)
      return
      end
