      program testvheo2
      implicit none
c
      real*8 R, z, theta, r1, v, vrccsdt, pi, abserr,
     +       relerr1, relerr2, tmp, vOO
      integer n, i
c
      pi=dacos(-1d0)
      abserr=0d0
      relerr1=0d0
      relerr2=0d0
      n=0
      write(*,*)'Using atomic units, theta in degrees'
      write(*,*)'Input format:  nr  R  theta   rOO   vHeO2'
      write(*,*)'The potential is printed when "nr" equals zero or'
      write(*,*)'a multiple of 150. "vHeO2" may be set to zero.'
10    continue
         read(*,*,end=100)i,R,theta,r1,vrccsdt
         n=n+1
         z=cos(theta*pi/180d0)
         call vheo2(R,z,r1,v)
         if (vrccsdt.lt.0d0) then
            tmp=dabs(v-vrccsdt)
            if (tmp.gt.abserr) abserr=tmp
         endif
         if (vrccsdt.gt.7.65d-4) then
            tmp=dabs(v/vrccsdt-1d0)
            if (tmp.gt.relerr1) relerr1=tmp
         endif
         if (R.gt.7d0) then
            tmp=dabs(v/vrccsdt-1d0)
            if (tmp.gt.relerr2) relerr2=tmp
         endif
         if (i.eq.150*(i/150)) then
            call xsigma(r1,vOO)
            write(*,'(1a,i3,3(1a,f10.6))')'nr=',i,' R=',R,' theta=',
     +      theta,' rOO=',r1
            write(*,200)'vHeO2 = ',vrccsdt,' (ab initio)'
            write(*,200)'vHeO2 = ',v,' (fit)'
            write(*,200)'vO2   = ',voo
            write(*,*)
200         format(1a,f18.12,1a)
         endif
         go to 10
100   continue
      write(*,*)'number of points ',n
      write(*,*)'Largest relative error for V>7.65e-4: ',relerr1*100d0,
     +' %'
      write(*,*)'Largest relative error for R>7:       ',relerr2*100d0,
     +' %'
      write(*,*)'Largest absolute error for V<0:       ',abserr
      end
c
************************************************************************
c
      subroutine vheo2(R,z,r1,v)
      implicit none
      real*8 R, z, r1, v
c
c input:
c     R: distance from He to center of mass of O2 (bohr).
c     z=cos(theta), theta=0 corresponds to the linear geometry.
c     r1: O-O bond distance (bohr).
c
c output:
c     v:  He + O2 interaction potential (hartree).
c
c This P.E.S. is based on RCCSD(T) calculations with the MOLPRO package,
c P. J. Knowles, C. Hampel and H.-J. Werner, J. Chem. Phys. 99 (1993) 5219.
c
c The grid used includes 750 points for which v<0.1 hartree in the region
c 3.05 <= R <= 15, 1.9 < r1 < 3.1
c For r1>4 and R<2.5 the potential may become unphysical
c The c_n,l with n=6, 8 are fixed from a fit to points with R>10.
c
c This subroutine will automatically perform a self-test on the first call,
c by checking the potential for one specific geometry.
c
c Gerrit C. Groenenboom, 21 August 1999.
c
      logical ifirst/.true./
      real*8 beta, betar, c(118)
      data beta/2.15d0/, betar/0.055d0/
      data c/
     +     4.09722719689086d+00,    -5.87850086171064d+00,
     +     2.62860977876725d+01,     5.33977970603514d+03,
     +    -7.03836915442398d+03,     3.08701418724401d+03,
     +    -4.49162257820320d+02,    -1.74812523863630d+04,
     +     2.29039582871963d+04,    -9.93850262047042d+03,
     +     1.42335345746197d+03,     1.77736241396470d+04,
     +    -2.33288030910082d+04,     1.01352371811363d+04,
     +    -1.45652788815395d+03,    -6.23418699018776d+03,
     +     8.01027388055036d+03,    -3.40685290682507d+03,
     +     4.78878899906915d+02,    -3.09385962612368d+03,
     +     4.11820103176701d+03,    -1.82194380085292d+03,
     +     2.69674770177646d+02,     1.12016768336068d+04,
     +    -1.46813538957560d+04,     6.38150949160146d+03,
     +    -9.16271316952770d+02,    -1.13367913344927d+04,
     +     1.48805598612998d+04,    -6.47556198430401d+03,
     +     9.31287308474901d+02,     3.96349644539987d+03,
     +    -5.09619063647082d+03,     2.17171045556329d+03,
     +    -3.06482889284378d+02,     5.36526094399032d+02,
     +    -7.10309442465364d+02,     3.14112121326586d+02,
     +    -4.66097444499953d+01,    -2.33873817730045d+03,
     +     3.05248585569763d+03,    -1.32364829295424d+03,
     +     1.89906454609946d+02,     2.37889179406219d+03,
     +    -3.11439979610650d+03,     1.35423641791217d+03,
     +    -1.95078301958074d+02,    -7.86519888471896d+02,
     +     1.00565170713353d+03,    -4.27199326425664d+02,
     +     6.02789202579294d+01,    -1.84135627390178d+01,
     +     2.42246635614066d+01,    -1.11610909492054d+01,
     +     1.77480604349836d+00,     1.45465852414389d+02,
     +    -1.87253678047965d+02,     8.10891496006745d+01,
     +    -1.17030768297363d+01,    -1.45733059813131d+02,
     +     1.89161213069541d+02,    -8.24141135939573d+01,
     +     1.19684422084845d+01,     4.59828585364982d+01,
     +    -5.79512931436981d+01,     2.44608084322934d+01,
     +    -3.45233993959894d+00,    -1.22239254764827d+01,
     +     6.99567706600854d+00,    -4.94111744139816d-01,
     +    -4.49529381223018d+00,     4.39827732944237d+00,
     +     6.69408361918515d-01,    -3.90839419951740d+02,
     +    -3.57488087850102d+02,     3.22042997461517d+02,
     +    -1.39598770751104d+02,    -4.18882366674752d+02,
     +     1.65799083233748d+02,     1.23172469401706d+01,
     +     1.46573401817661d+02,    -6.64328930023274d+01,
     +    -2.84334402957489d+04,     2.25544569447039d+05,
     +    -1.93364408972475d+05,     5.12862150223910d+04,
     +    -7.51711206190477d+04,     3.64392258273544d+05,
     +    -3.03121711252478d+05,     8.73346907324946d+04,
     +    -3.45519750882223d+04,     4.25719100787138d+05,
     +    -3.75943308872289d+05,     9.45902804164792d+04,
     +    -4.60944149675849d+03,     1.14916075693361d+05,
     +    -1.04417118162543d+05,     2.54022051744389d+04,
     +     2.88402408686071d+06,    -2.71785831822367d+07,
     +     2.52655446346514d+07,    -7.03236130965375d+06,
     +     3.11298177116258d+06,    -3.63253073945702d+07,
     +     3.45445576187005d+07,    -9.27396845434792d+06,
     +     1.74021246580691d+06,    -4.13399996654562d+07,
     +     3.66427868291820d+07,    -8.60763269788007d+06,
     +    -8.22733314352684d+04,    -6.46296378998269d+06,
     +     5.79777234997689d+06,    -1.26247393590976d+06,
     +    -6.29973251125557d+04,     2.17132361189737d+05,
     +    -1.87680229897387d+05,     6.28352247097358d+04/
c
      real*8 Ra, za, Rb, zb, pla(4), plb(4), pl(9),
     + Ran(4), Rbn(4), r1n(4), sum, Rn(4), f_long,
     + Rsave, zsave, r1save, v0
      integer i, j, k, ind, kmax
c
      if (ifirst) then
         write(*,'(1a)')
     +'=================================================',
     +'RCCSD(T) He+O2 interaction potential,',
     +'Gerrit C. Groenenboom and Izabela M. Struniewicz,',
     +'J. Chem. Phys. vol 113, December (2000)',
     +'E-mail: gerritg@theochem.kun.nl',
     +'================================================='
c
c       Save the point
c
        Rsave=R
        zsave=z
        r1save=r1
c
c       ... and first check whether this point gives the right value:
c
        R=7d0
        z=0.5d0
        r1=2.5d0
        v0=-0.866042743602981517D-04
        go to 20
      endif
10    continue
      if (ifirst) then
        if (dabs((v/v0)-1d0).gt.1d-12) then
           write(*,*)'** ERROR in vHeO2: integrity check not passed'
           write(*,*)'** v(',R,',',z,',',r1,')=',v
           write(*,*)'** this should be ',v0
           stop 230
        endif
c       OK, check passed, do the required point
        R=Rsave
        z=zsave
        r1=r1save
        ifirst=.false.
      endif
c
20    continue
      v=0d0
c
      if (dabs(z).gt.1d0) then
         write(*,*)'** ERROR in VHEO2: illegal cos(theta) = z = ',z
         stop 231
      endif
      call jac2yu(R,z,r1,Ra,za,Rb,zb)
      call p_leg(pla,3,za)
      call p_leg(plb,3,zb)
      Ran(1)=dexp(-beta*Ra)
      Rbn(1)=dexp(-beta*Rb)
      r1n(1)=dexp(-betar*r1**3)
      do i=2,4
        Ran(i)=Ra*Ran(i-1)
        Rbn(i)=Rb*Rbn(i-1)
        r1n(i)=r1*r1n(i-1)
      end do
c
c     the r1 independent part of the short range
c
      ind=0
      do i=1,3
         ind=ind+1
         v=v+(Ran(1)*pla(i)+Rbn(1)*plb(i))*c(ind)
      end do
c
c     the r1 dependent part of the short range
c
      do i=1,4
         do j=1,4
            sum=0d0
            do k=1,4
               ind=ind+1
               sum=sum+c(ind)*r1n(k)
            end do
            v=v+(Ran(i)*pla(j)+Rbn(i)*plb(j))*sum
         end do
      end do
c
c     the long range
c
      do i=1,4
         Rn(i)=f_long(R,beta,4+2*i)
      end do
      call p_leg(pl,8,z)
c
      kmax=2
      do i=1,4
         if (i.eq.3) kmax=kmax+1
         do j=1,i+1
            ind=ind+1
            sum=c(ind)
            do k=1,kmax
               ind=ind+1
               sum=sum+r1n(k)*c(ind)
            end do
            v=v+Rn(i)*pl(2*j-1)*sum
         end do
      end do

      if (ifirst) go to 10
      end
c
************************************************************************
c
      subroutine jac2yu(R,z,r1,Ra,za,Rb,zb)
      implicit none
      real*8 R, r1, z, Ra, za, Rb, zb
      real*8 alpha, tmpa, tmpb, x, y2

      if (R.eq.0d0) then
         Ra=0.5d0*r1
         za=-1d0
         Rb=Ra
         zb=za
      else
         alpha=r1/(2d0*R)
         tmpa=dsqrt(1d0-2d0*alpha*z+alpha*alpha)
         tmpb=dsqrt(1d0+2d0*alpha*z+alpha*alpha)
         if (tmpa.eq.0d0) then
            za=1d0
         else
            za=(z-alpha)/tmpa
         end if
         if (tmpb.eq.0d0) then
            zb=1d0
         else
            zb=-(z+alpha)/tmpb
         end if
         x=R*z
         y2=R*R*(1d0-z*z)
         Ra=dsqrt( ( x-0.5d0*r1)**2 + y2 )
         Rb=dsqrt( ( x+0.5d0*r1)**2 + y2 )
      endif
      end
c
************************************************************************
c
      real*8 function f_long(r,beta,n)
      implicit none
      real*8 r, beta
      integer n
c
c     f_long(r,beta,n)=d_n(beta*r,n)/r^n
c
c     Tang-Toennies damping function:
c     d_n(x,n) = 1 - exp(-x) * sum_k=0:n x^k / k!
c
c     For small x (when d_n<1e-8) use:
c     d_n(x,n) = exp(-x)* sum_k=n+1:infty x^k/k!
c
      real*8  sum, term, dex, rn, dn, x
      integer i
c
      if (r.eq.0d0) then
         f_long=0d0
         return
      endif
c
      x=beta*r
      sum=1d0
      term=1d0
      rn=1d0
      do i=1,n
         rn=rn*r
         term=term*x/dfloat(i)
         sum=sum+term
      end do
      dex=dexp(-x)
      dn=1d0-dex*sum
      if (dabs(dn).lt.1d-8) then
c
c        alternative formula for small x
c
         sum=0d0
         do i=n+1,1000
            term=term*x/dfloat(i)
            sum=sum+term
            if (term/sum .lt. 1d-8) go to 10
         enddo
         write(*,*)'** ERROR in F_LONG: no convergence, x=',x,' n=',n
         stop 232
10       continue
         dn=sum*dex
      endif
      f_long=dn/rn
      end
c
************************************************************************
c
      subroutine p_leg(p,lmax,z)
      implicit none
      real*8 p(0:lmax), z, di
      integer lmax, i
c
      p(0) = 1.d0
      if (lmax.eq.0) return
      p(1) = z
      di = 0d0
      do i=1,lmax-1
        di = di + 1d0
        p(i+1) = ( (2.d0*di + 1.d0) * z * p(i)
     +                - di * p(i-1) ) / (di + 1d0)
      end do
      end
c
************************************************************************
c
      subroutine xsigma(r,v)
      parameter ( npts = 24000 )
      implicit double precision(a-h,o-z)
c     r & v in atomic units.
c
c  ground state x3sgm potential energy curve of 02
c  rkr values from cosby and helms for 1.79480 < r < 4.0840
c  dissociation energy 41260.2 cm-1 from brix and herzberg
c  zero point energy 788.1 cm-1 from cosby and helms
c  short range from ron friedman
c  long range  -17.57/r**6 assumed valid for r > 5.0
c  nine points 4.2< r < 5.0 added to join smoothly to long range
c  generated by arthur allison july 1991 using difference between
c  4.5 and 5.0 from steve guberman
c
      dimension ra(200),pr(200),c(4,200)
      data nc,yp0,ypm/148,-1.018436615d0,1.3494d-03/
      data (ra(i),i=1,76)/
     *  0.17948246d+01, 0.17963175d+01, 0.17978482d+01, 0.17993977d+01,
     *  0.18010229d+01, 0.18027425d+01, 0.18045189d+01, 0.18064086d+01,
     *  0.18083550d+01, 0.18103959d+01, 0.18125313d+01, 0.18147234d+01,
     *  0.18169911d+01, 0.18193154d+01, 0.18217154d+01, 0.18241720d+01,
     *  0.18267043d+01, 0.18292743d+01, 0.18319199d+01, 0.18346411d+01,
     *  0.18374001d+01, 0.18402347d+01, 0.18431449d+01, 0.18461118d+01,
     *  0.18491731d+01, 0.18522912d+01, 0.18554848d+01, 0.18587729d+01,
     *  0.18621556d+01, 0.18656138d+01, 0.18691664d+01, 0.18728136d+01,
     *  0.18765742d+01, 0.18804292d+01, 0.18843976d+01, 0.18884794d+01,
     *  0.18926746d+01, 0.18969832d+01, 0.19014241d+01, 0.19060161d+01,
     *  0.19107404d+01, 0.19156159d+01, 0.19206426d+01, 0.19258393d+01,
     *  0.19312251d+01, 0.19367998d+01, 0.19425823d+01, 0.19485916d+01,
     *  0.19548655d+01, 0.19613851d+01, 0.19682070d+01, 0.19753502d+01,
     *  0.19828524d+01, 0.19907515d+01, 0.19990851d+01, 0.20079102d+01,
     *  0.20173021d+01, 0.20273177d+01, 0.20380891d+01, 0.20497109d+01,
     *  0.20623721d+01, 0.20763183d+01, 0.20918518d+01, 0.21095208d+01,
     *  0.21301944d+01, 0.21555923d+01, 0.21902310d+01, 0.22157990d+01,
     *  0.22816991d+01, 0.23539947d+01, 0.23859878d+01, 0.24333632d+01,
     *  0.24715168d+01, 0.25049461d+01, 0.25353896d+01, 0.25637355d+01/
      data (ra(i),i=77,148)/
     *  0.25905507d+01, 0.26161565d+01, 0.26408174d+01, 0.26647036d+01,
     *  0.26879472d+01, 0.27106617d+01, 0.27329227d+01, 0.27547868d+01,
     *  0.27763297d+01, 0.27975891d+01, 0.28186029d+01, 0.28393899d+01,
     *  0.28599879d+01, 0.28804348d+01, 0.29007493d+01, 0.29209316d+01,
     *  0.29410005d+01, 0.29609938d+01, 0.29809115d+01, 0.30007726d+01,
     *  0.30205769d+01, 0.30403623d+01, 0.30601100d+01, 0.30798765d+01,
     *  0.30996431d+01, 0.31194285d+01, 0.31392328d+01, 0.31590939d+01,
     *  0.31790116d+01, 0.31990049d+01, 0.32190927d+01, 0.32392750d+01,
     *  0.32595517d+01, 0.32799797d+01, 0.33005399d+01, 0.33212513d+01,
     *  0.33421328d+01, 0.33632032d+01, 0.33845005d+01, 0.34060055d+01,
     *  0.34277752d+01, 0.34497905d+01, 0.34721271d+01, 0.34947660d+01,
     *  0.35177640d+01, 0.35411588d+01, 0.35649882d+01, 0.35892901d+01,
     *  0.36141211d+01, 0.36395569d+01, 0.36656540d+01, 0.36924881d+01,
     *  0.37201915d+01, 0.37488587d+01, 0.37786218d+01, 0.38096323d+01,
     *  0.38420978d+01, 0.38762451d+01, 0.39123200d+01, 0.39506815d+01,
     *  0.39917074d+01, 0.40359270d+01, 0.40839828d+01, 0.42d+01,
     *  0.43d+01,       0.44d+01,       0.45d+01,       0.46d+01,
     *  0.47d+01,       0.48d+01,       0.49d+01,       0.50d+01/
      data (pr(i),i=1,76)/
     * -0.12174598d-01,-0.13688451d-01,-0.15266193d-01,-0.16903654d-01,
     * -0.18597132d-01,-0.20343372d-01,-0.22139485d-01,-0.23982944d-01,
     * -0.25871496d-01,-0.27803170d-01,-0.29776206d-01,-0.31789060d-01,
     * -0.33840347d-01,-0.35928850d-01,-0.38053471d-01,-0.40213235d-01,
     * -0.42407253d-01,-0.44634729d-01,-0.46894947d-01,-0.49187246d-01,
     * -0.51511020d-01,-0.53865722d-01,-0.56250838d-01,-0.58665903d-01,
     * -0.61110484d-01,-0.63584175d-01,-0.66086604d-01,-0.68617436d-01,
     * -0.71176354d-01,-0.73763069d-01,-0.76377315d-01,-0.79018852d-01,
     * -0.81687464d-01,-0.84382962d-01,-0.87105179d-01,-0.89853955d-01,
     * -0.92629171d-01,-0.95430714d-01,-0.98258492d-01,-0.10111243d+00,
     * -0.10399248d+00,-0.10689857d+00,-0.10983069d+00,-0.11278879d+00,
     * -0.11577288d+00,-0.11878293d+00,-0.12181893d+00,-0.12488087d+00,
     * -0.12796875d+00,-0.13108257d+00,-0.13422229d+00,-0.13738792d+00,
     * -0.14057943d+00,-0.14379683d+00,-0.14704006d+00,-0.15030915d+00,
     * -0.15360407d+00,-0.15692484d+00,-0.16027149d+00,-0.16364407d+00,
     * -0.16704269d+00,-0.17046749d+00,-0.17391871d+00,-0.17739663d+00,
     * -0.18090168d+00,-0.18443440d+00,-0.18799547d+00,-0.18978691d+00,
     * -0.19158624d+00,-0.18978691d+00,-0.18799547d+00,-0.18443440d+00,
     * -0.18090168d+00,-0.17739663d+00,-0.17391871d+00,-0.17046749d+00/
      data (pr(i),i=77,148)/
     * -0.16704269d+00,-0.16364407d+00,-0.16027149d+00,-0.15692484d+00,
     * -0.15360407d+00,-0.15030915d+00,-0.14704006d+00,-0.14379683d+00,
     * -0.14057943d+00,-0.13738792d+00,-0.13422229d+00,-0.13108257d+00,
     * -0.12796875d+00,-0.12488087d+00,-0.12181893d+00,-0.11878293d+00,
     * -0.11577288d+00,-0.11278879d+00,-0.10983069d+00,-0.10689857d+00,
     * -0.10399248d+00,-0.10111243d+00,-0.98258492d-01,-0.95430714d-01,
     * -0.92629171d-01,-0.89853955d-01,-0.87105179d-01,-0.84382962d-01,
     * -0.81687464d-01,-0.79018852d-01,-0.76377315d-01,-0.73763069d-01,
     * -0.71176354d-01,-0.68617436d-01,-0.66086604d-01,-0.63584175d-01,
     * -0.61110484d-01,-0.58665903d-01,-0.56250838d-01,-0.53865722d-01,
     * -0.51511020d-01,-0.49187246d-01,-0.46894947d-01,-0.44634729d-01,
     * -0.42407253d-01,-0.40213235d-01,-0.38053471d-01,-0.35928850d-01,
     * -0.33840347d-01,-0.31789060d-01,-0.29776206d-01,-0.27803170d-01,
     * -0.25871496d-01,-0.23982944d-01,-0.22139485d-01,-0.20343372d-01,
     * -0.18597132d-01,-0.16903654d-01,-0.15266193d-01,-0.13688451d-01,
     * -0.12174598d-01,-0.10729373d-01,-0.93581104d-02,-0.690d-02,
     * -0.500d-02,     -0.357d-02,     -0.260d-02,     -0.207d-02,
     * -0.170d-02,     -0.145d-02,     -0.127d-02,     -0.11245d-02/
      logical first/.true./
      if (first) then
         write(*,'(1a)')
     +'=========================================',
     +'O2 potential, J. F. Babb and A. Dalgarno,',
     +'Phys. Rev. A, 51, 3021 (1995)',
     +'========================================='
         first=.false.
      endif
c
c     asym used only for short range fit a*exp(-alpha*r) rel v11 = 0
      asym = 41260.2d0/219474.4354d0
      nint = 1
      call splicn(ra,pr,nc,c,yp0,ypm)
c   x state of o2
        if ( r .gt. ra(1) ) then
          if ( r .lt. ra(nc) ) then
c           spline fit
   70       if ( (r-ra(nint+1)) .lt. 0 )  goto 71
              nint = nint + 1
              goto 70
   71       continue
            rl = r - ra(nint)
            ru = ra(nint+1) - r
            v = c(4,nint)*rl    + c(3,nint)*ru
     *           + c(2,nint)*rl**3 + c(1,nint)*ru**3
          else
c           long range
            v = -17.57d0/r**6
          endif
        else
c         short range
          v = 5757.005316d0*dexp(-5.792467035d0*r) - asym
        endif
c       write(21,*)r,v
c       ucm = 42048.3 + v*219474.4354
      return
      end
c
      subroutine splicn(x,y,m,c,ypo,ypm)
      implicit double precision(a-h,o-z)
c
c     spline function subroutine
c     set of m pts,abscissae in x,ordinates in y,initial deriv yp0,
c     final deriv ypm, coeffs left in c  max of 200 points
c
      dimension x(m),y(m),d(200),p(200),e(200),c(4,m),a(200,3),
     *       b(200),z(200)
c
      mm = m-1
      do 2 k = 1,mm
        d(k) = x(k+1) - x(k)
        p(k) = d(k)/6.0d0
        e(k) = (y(k+1) - y(k))/d(k)
    2 continue
      p(m) = 0.0d0
      b(1) = (e(1) - ypo)/(2.0d0*p(1))
      do 3 k = 2,mm
        b(k) = e(k) - e(k-1)
    3 continue
c
      b(m) = ypm - e(m-1)
      a(1,3) = 0.5d0
      do 4 k = 2,m
        a(k,2) = 2.0d0*(p(k-1) + p(k)) - p(k-1)*a(k-1,3)
        b(k)   = b(k) - p(k-1)*b(k-1)
        a(k,3) = p(k)/a(k,2)
        b(k)   = b(k)/a(k,2)
    4 continue
c
      z(m) = b(m)
      do 6 i = 1,mm
        k    = m-i
        z(k) = b(k) - a(k,3)*z(k+1)
    6 continue
c
      do 7 k = 1,mm
        q   = 1.0d0/(6.0d0*d(k))
        c(1,k) = z(k)*q
        c(2,k) = z(k+1)*q
        c(3,k) = y(k)/d(k) - z(k)*p(k)
        c(4,k) = y(k+1)/d(k) - z(k+1)*p(k)
    7 continue
      return
      end

