cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc
      subroutine LiHF_BO_init()
      call prepot_LiHF_BO
      return
      end
cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc
      subroutine LiHF_BO_pes(r,v)
      implicit real*8 (a-h,o-z)
      dimension r(3)
      data auang,aucclm,vbase/.52917706d0,1.59359d-3,141.2d0/
      rab=auang*r(1)
      rbc=auang*r(3)
      rac=auang*r(2)
      call potbo_LiHF_BO(rab,rbc,rac,v)
      v=(v+141.2d0)*aucclm
      return
      end
c----------------------------------------------------------------------c
      subroutine potbo_LiHF_BO(rab,rbc,rac,v)
c
      implicit real*8(a-h,o-z)
      common/b21_LiHF_BO/ de1,de2,de3,re1,re2,re3
     >   ,b1,b2,b3
       common/b22_LiHF_BO/a11,a12,a13,a21,a22,a23,a31,a32,a33,a14,a24,
     $a34
       common/btanh_LiHF_BO/ref1,ref2,ref3,sl1,sl2,sl3,r01,r02,r03
       common/repuls_LiHF_BO/a,n1,n2,n3
       common/refine_LiHF_BO/agb,cgbx,cgby,cgbz,agw,cgwx,cgwy,cgwz,
     >               xn0b,yn0b,zn0b,xn0w,yn0w,zn0w
c
      x=rab
      y=rbc
      z=rac
c
      if(x.le.0.d0)x=1.d-4
      if(y.le.0.d0)y=1.d-4
      if(z.le.0.d0)z=1.d-4
c
c hyperbolic tangent damping functions
c
c      w1=.5d0*(1d0+tanh(sl1*(x-ref1)))
c      w2=.5d0*(1d0+tanh(sl2*(y-ref2)))
c      w3=.5d0*(1d0+tanh(sl2*(z-ref3)))
c      ww=w1*w2*w3
c
c  diatomic terms and their corrections
c
      v=0.d0
       ex=exp(-b1*(x-re1))
       ey=exp(-b2*(y-re2))
       ez=exp(-b3*(z-re3))
       vxx=-de1*ex*(a11+ex*(a12+ex*(a13+ex*a14)))
       vyy=-de2*ey*(a21+ey*(a22+ey*(a23+ey*a24)))
       vzz=-de3*ez*(a31+ez*(a32+ez*(a33+ez*a34)))
       tx=de1*ex**n1
       ty=de2*ey**n2
       tz=de3*ez**n3
       vxx=vxx+a*tx
       vyy=vyy+a*ty
       vzz=vzz+a*tz
       vbiat=vxx+vyy+vzz
c       write(6,*)n1,n2,n3,tx,ty,tx
c
c  the triatomic term and its corrections
c
        call pol3_LiHF_BO(x,y,z,v3xyz)
c
c      if(x.lt.re1)x=re1
c      if(y.lt.re2)y=re2
c      if(z.lt.re3)z=re3
c      call pol3(x,y,z,v3re)
c
c      vtriat=ww*v3xyz+(1.d0-ww)*v3re
       vcbar=exp(-cgbx*(ex-xn0b)*(ex-xn0b))*
     > exp(-cgby*(ey-yn0b)*(ey-yn0b))*
     > exp(-cgbz*(ez-zn0b)*(ez-zn0b))
       vcwel=exp(-cgwx*(ex-xn0w)*(ex-xn0w))*
     > exp(-cgwy*(ey-yn0w)*(ey-yn0w))*
     > exp(-cgwz*(ez-zn0w)*(ez-zn0w))
       vtriat=.9935d0*v3xyz + agb*vcbar + agw*vcwel
c      if(vcbar.gt..001d0.or.vcwel.gt..001d0)write(6,*)agb,vcbar,agw,vcwel
       v=vbiat+vtriat
       if(v.gt.331.d0)v=331.0d0+log(v-330.0d0)
       vshow=v+de2
c     write(8,11)x,y,z,sx,sy,sz,vxx,vyy,vzz,vtriat,vshow
11      format(6f5.2,5f8.2)
                             !/ NO RETURN in original code
      RETURN
      end
c----------------------------------------------------------------------c
      subroutine pol3_LiHF_BO(x,y,z,v3)
      implicit real*8(a-h,o-z)
c
       common/b21_LiHF_BO/de1,de2,de3,re1,re2,re3,b1,b2,b3
      common/b23_LiHF_BO/coef(65),iexp(65),jexp(65),kexp(65),nt
      dimension tx(6),ty(6),tz(6)
c
      tx(1)=1.d0
      ty(1)=1.d0
      tz(1)=1.d0
       ex=exp(-b1*(x-re1))
       ey=exp(-b2*(y-re2))
       ez=exp(-b3*(z-re3))
      do 10 i=2,6
      ii=i-1
       tx(i)=tx(ii)*ex
       ty(i)=ty(ii)*ey
       tz(i)=tz(ii)*ez
   10 continue
      w=0.d0
      do 20 i=1,nt
      ix=iexp(i)
      jx=jexp(i)
      kx=kexp(i)
      w=w+coef(i)*tx(ix)*ty(jx)*tz(kx)
   20 continue
      v3=w
      return
      end
c----------------------------------------------------------------------c
*deck potbo
c
c  this deck potbo contains:
c            the interface routine pot
c                       vbase adjusts the energy zero
c                       convr converts internuclear distances to a
c                       convv converts potential values to kcal.mol
c                       r are internucl. distances in the order rac,rab,rbc
c                       vv is the output potential value in desired units
c            the routine prepot for reading potential parameters
c            the routine potbo for calculating bo potentials
c               as a function of rab,rbc,rac (in a). it returns as v the
c               value of the potential in kcalmol. the energy zero is set
c               at total dissociation.
c           the soubroutine pol3 for calculating the bo polynomial
c
c
      subroutine prepot_LiHF_BO
c
      implicit real*8(a-h,o-z)
      common/b21_LiHF_BO/ de1,de2,de3,re1,re2,re3
     >   ,b1,b2,b3
      common/b23_LiHF_BO/cof(65),ie(65),je(65),ke(65),nt
      dimension coef(65),iexp(65),jexp(65),kexp(65)
      common/b22_LiHF_BO/a11,a12,a13,a21,a22,a23,a31,a32,a33
     >   ,a14,a24,a34
      common/b25_LiHF_BO/ pg
       common/btanh_LiHF_BO/ref1,ref2,ref3,sl1,sl2,sl3,r01,r02,r03
       common/repuls_LiHF_BO/a,n1,n2,n3
       common/refine_LiHF_BO/agb,cgbx,cgby,cgbz,agw,cgwx,cgwy,cgwz,
     >               xn0b,yn0b,zn0b,xn0w,yn0w,zn0w
c
       data coef /
     *   0.9224052d+03,  0.3272860d+03,  0.3381706d+03, -0.1048643d+04,
     *  -0.2962569d+03, -0.3839140d+03, -0.7478068d+03, -0.7533124d+03,
     *  -0.1715470d+03, -0.2134407d+01,  0.5508481d+03,  0.1134382d+04,
     *  -0.1160042d+04,  0.3557390d+03,  0.1017086d+02,  0.6969160d+03,
     *   0.1258877d+04,  0.5197582d+02,  0.4328018d+03, -0.1062101d+03,
     *   0.5471942d+02, -0.3211568d+02, -0.1333515d+03, -0.4993582d+03,
     *  -0.3928734d+03,  0.1143333d+04, -0.2252141d+03,  0.2641766d+03,
     *   0.9964729d+02, -0.3031654d+03, -0.9919408d+03,  0.7123589d+02,
     *  -0.5026490d+03, -0.2019268d+03,  0.1085769d+03, -0.3689024d+02,
     *   0.1259192d+03,  0.1398510d+02, -0.5489551d+01, -0.2721741d+02,
     *  -0.8672702d+01,  0.9243725d+02,  0.9455495d+02,  0.4378558d+02,
     *  -0.3012895d+03,  0.5875433d+02, -0.7644942d+01, -0.1544242d+03,
     *   0.1174672d+01,  0.1881099d+02,  0.2699340d+03, -0.3803897d+02,
     *   0.2211474d+02,  0.1713199d+03,  0.6709845d+02, -0.2999395d+02,
     *   0.1285414d+02,  0.5457174d+02, -0.9745477d+02, -0.7197123d+01,
     *  -0.1377454d+01, -0.2120028d+02,  0.2071087d+02,  0.7003032d+00,
     *   0.4229195d+01/
       data iexp/
     * 1, 1, 0, 2, 1, 2, 1, 0, 1, 0, 3, 2, 1, 3, 2, 1, 0, 2, 1, 0,
     * 1, 0, 4, 3, 2, 1, 4, 3, 2, 1, 0, 3, 2, 1, 0, 2, 1, 0, 1, 0,
     * 5, 4, 3, 2, 1, 5, 4, 3, 2, 1, 0, 4, 3, 2, 1, 0, 3, 2, 1, 0,
     * 2, 1, 0, 1, 0/
       data jexp/
     * 1, 0, 1, 1, 2, 0, 1, 2, 0, 1, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2,
     * 0, 1, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 0, 1, 2, 0, 1,
     * 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 2, 3,
     * 0, 1, 2, 0, 1/
       data kexp/
     * 0, 1, 1, 0, 0, 1, 1, 1, 2, 2, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2,
     * 3, 3, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4,
     * 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3,
     * 4, 4, 4, 5, 5/
c
       data vvre,vvr0,eps/1d-4,.01d0,.001d0/
       data tailbx,tailby,tailbz/.111d0,.037d0,.078d0/
       data tailwx,tailwy,tailwz/.088d0,.096d0,.042d0/
       data x0b,y0b,angb/1.6523d0,1.3715d0,72.5235d0/
       data x0w,y0w,angw/1.8848d0,0.9363d0,106.929d0/
       agw = -1.8d0
       agb = -0.1d0
c
    2 format(i2)
    3 format(e15.7,3i3)
c
      pg=3.1415926536d0
c
c   dissociation energies
       de1=137.59d0
       de2=141.20d0
       de3=58.d0
c
c   b parameters
c
       b1=.970755d0
       b2=2.194185d0
       b3=1.17086d0
c
c   equilibrium distances
      re1=1.563864d0
      re2=0.916808d0
      re3=1.5957d0
c
c    diatomic parameters
       a11=2.30439d0
       a12=-2.30172d0
       a13=1.69028d0
       a14=-.69295d0
       a21=2.07810d0
       a22=-1.25667d0
       a23=0.27905d0
       a24=-0.10047d0
       a31=2.15123d0
       a32=-1.38219d0
       a33=0.31067d0
       a34=-.07972d0
c
c     finds the largest short range distance at which diatomic
c     potentials are zero (r01,r02,r03)
      do 1 i=1,3
       if(i.ne.1)goto 63
       b=b1
       re=re1
       de=de1
       a1=a11
       a2=a12
       a3=a13
       a4=a14
       goto67
63     continue
       if(i.ne.2)goto65
       b=b2
       re=re2
       de=de2
       a1=a21
        a2=a22
       a3=a2e
        a4=a24
        goto 67
65      continue
       b=b3
       re=re3
       de=de3
       a1=a31
       a2=a32
       a3=a33
       a4=a34
67     continue
       dx=-.1d0
       x=re-dx
70     continue
       x=x+dx
       ex=exp(-b*(x-re))
       v=-de*ex*(a1+ex*(a2+ex*(a3+ex*a4)))
       if(v.lt.0d0)goto70
       if(abs(dx).lt.0.0001d0)goto110
        x=x-dx
        dx=.1d0*dx
       goto70
110    continue
       if(i.eq.1)r01=x
       if(i.eq.2)r02=x
       if(i.eq.3)r03=x
1      continue
c
c    determines the value at which tanh has a flex
       ref1=.49d0*r01+.51d0*re1
       ref2=.49d0*r02+.51d0*re2
       ref3=.49d0*r03+.51d0*re3
c
c   determines the slope of tanh at the flex
       alneps=log(1d0/eps-1d0)
       sl1=.5d0*alneps/(re1-ref1)
       sl2=.5d0*alneps/(re2-ref2)
       sl3=.5d0*alneps/(re3-ref3)
c
c  determines the parameters for the additional repulsive inner
c  diatomic term the coefficient a and the powers n1,n2,n3
c
       a=vvre
       n1=log(vvr0/vvre)/(-b1*(r01-re1))
       n2=log(vvr0/vvre)/(-b2*(r02-re2))
       n3=log(vvr0/vvre)/(-b3*(r03-re3))
c
c  determines the parameters for the additional repulsive short range
c  terms of diatoms (given as D*a*BO**n)
c
       a=vvre
       n1=log(vvr0/vvre)/(-b1*(r01-re1))
       n2=log(vvr0/vvre)/(-b2*(r02-re2))
       n3=log(vvr0/vvre)/(-b3*(r03-re3))
c
c  calculates the parameters for gaussian bond order corrections to
c  the barrier and the well
c
       gamb=angb*pg/180.d0
       gamw=angw*pg/180.d0
       z0b=sqrt(x0b*x0b+y0b*y0b-2.d0*x0b*y0b*cos(gamb))
       z0w=sqrt(x0w*x0w+y0w*y0w-2.d0*x0w*y0w*cos(gamw))
       xn0b=exp(-b1*(x0b-re1))
       yn0b=exp(-b2*(y0b-re2))
       zn0b=exp(-b3*(z0b-re3))
       xn0w=exp(-b1*(x0w-re1))
       yn0w=exp(-b2*(y0w-re2))
       zn0w=exp(-b3*(z0w-re3))
       cl1=log(.5d0)
       cgbx=-cl1/(tailbx*tailbx)
       cl1=log(.5d0)
       cgby=-cl1/(tailby*tailby)
       cl1=log(.5d0)
       cgbz=-cl1/(tailbz*tailbz)
       cl1=log(.5d0)
       cgwx=-cl1/(tailwx*tailwx)
       cl1=log(.5d0)
       cgwy=-cl1/(tailwy*tailwy)
       cl1=log(.5d0)
       cgwz=-cl1/(tailwz*tailwz)
c
c  coefficients for the three-body term
c
       nt=65
      do 10 i=1,nt
      cof(i)=coef(i)
      ie(i)=1+iexp(i)
      je(i)=1+jexp(i)
      ke(i)=1+kexp(i)
c     write(6,*)'prepot', cof(i),ie(i),je(i),ke(i)
   10 continue
c
      return
      end
