c
c  f77 -o vlam_h3.exe vlam_h3.f bkmp1.f bkmp2.f dmbe.f lsth2.f h3aguado.f
c  ./vlam_h3.exe
c
c  takes about 4 minutes on [Intel Pentium III 450 MHz, Linux, fort77]
c
      program vlam_h3
c     ---------------
c
      implicit double precision (a-h,o-z)
c
      parameter ( dvtry = 1.d-8, dvptry = 1.d-7 )
c
      parameter ( nmax = 65536, nstart = 8 )
c
      dimension vi(0:nmax,3)
c
      character*80 cselect
c
      character*100 cfile
c
      character*5 cend(6)
c
      character*1 conf
c
      parameter ( nri = 5 )
c
      dimension rhh(nri), vlam(nri,3), v(3), vp(3)
c
      dimension e(3), dVtot(3)
c
      COMMON /POTCM/ ri(3), POTE1, DH3L(3)
c
      LOGICAL DERIVL
      COMMON/DERCOM/DERIVL
c
c  Names specifying the six H3 surfaces:
c
      data cend / 'BKMP1', 'BKMP2', 'DMBE ', 'LSTH ',
     $     'Agua7', 'Agu11' /
c
c ***** The five default H2-molecule sizes:
c
      data rhh / 0.9d0, 1.28d0, 1.449d0, 1.618d0, 2.d0 /
c
      DERIVL = .FALSE.
c
      call prepot
c
      pi = acos(-1.d0)
      thninv = 1.d0 / nmax
c
cLegendre coefficents for BMKP2, BKMP, DMBE, LSTH, Agua7, & Agu11 H3 surfaces
cas a function of R, for H2-molecule sizes r in the range [0.6,6.0]
c
      write(6,2)
 2    format(/' Legendre coefficents for BMKP2, BKMP, DMBE,'
     $     ' LSTH, Agua7, & Agu11 H3 surfaces'/
     $     ' as a function of R, for H2-molecule sizes r'
     $     ' in the range [0.6,6.0]'//
     $     ' Give H2-size-indicator string to insert',
     $     ' in filenames (blank = default case):')
      read(5,'(a80)') cselect
c
      last = 80
      do while ( last .gt. 1 .and. cselect(last:last) .eq. ' ' )
         last = last - 1
      enddo
c
      if ( last .eq. 1 .and. cselect(last:last) .eq. ' ' ) then
c
         last = 0
         write(6,8) (rhh(i),i=1,5)
 8       format(/' r = {',5f6.3,' } (default): continue...')
c
      else
c
         do k = 1, last
            if ( cselect(k:k) .eq. ' ' ) stop
     $           ' STOP -- Error: embedded blanks in string. '
         enddo
c
         write(6,4)
 4       format(/' Give five H2-molecule sizes r: ',$)
         read(5,*) (rhh(i),i=1,5)
c
         do i = 1, 5
            if ( rhh(i) .lt. 0.6d0 .or. rhh(i) .gt. 6.d0 ) stop
     $           ' STOP -- Error: bad H2-molecule size. '
         enddo
c
         write(6,6) cselect(:last), (rhh(i),i=1,5)
 6       format(/' H2-size-indicator string = ',a/
     $        ' r = {',5f10.6,' } : Confirm? (y/n): ',$)
         read(5,'(a1)') conf
c
         if ( conf .ne. 'y' .and. conf .ne. 'Y' ) stop
     $        ' STOP -- Run cancelled. '
c
      endif
c
      do isurf = 6, 1, -1
c
         if ( last .eq. 0 ) write(6,
     $        '(/" Create",3(" v",i1,"_h3_",a,".txt"))')
     $        ( lam, cend(isurf)(:5-isurf/3+isurf/5+isurf/6),
     $        lam = 0, 4, 2 )
c
         do lam = 0, 4, 2
c
            k = max( 1 , min( 3 , lam ) )
c
c  Create the output filename for this lambda and surface:
c
            if ( last .eq. 0 ) then
               write(cfile,'("v",i1,"_h3_",a,".txt")') lam,
     $              cend(isurf)(:5-isurf/3+isurf/5+isurf/6)
            else
               write(cfile,'("v",i1,"_h3_",a,"_",a,".txt")')
     $              lam, cselect(:last),
     $              cend(isurf)(:5-isurf/3+isurf/5+isurf/6)
               if ( k .eq. 1 ) then
                  write(6,'(/" Create ",a)') cfile(:last+16)
               else
                  write(6,'("        ",a)') cfile(:last+16)
               endif
            endif
c
c  and open this file to write to it:
c
            open(k,file=cfile,form='FORMATTED',status='UNKNOWN')
c
            write(k,10) cend(isurf), lam, lam,
     $           ( lam, rhh(i), i = 1, nri )
 10         format('# ',a5,' lambda=',i1,': V_',i1,
     $           '(H3) [5 H2-molecule sizes]:',
     $           ' r, R in bohrs, V in hartrees.'/
     $           '#'/'#'/
     $           '#  R ',5(' V',i1,'(r=',f7.5,')')/
     $           '#==== ------------- -------------',
     $           ' ============= ------------- -------------')
c
         enddo
c
         jr = 16
c
         do while ( jr .le. 200 )
c
            R = 0.1d0 * jr
c
            nthmax = 0
            dm = 0.d0
c
            do ir = 1, nri
c
               rh = rhh(ir)
               rh_o_2 = 0.5d0 * rh
c
c  Avoid error halt in BKMP2 surface (and errors in other surfaces) at small r2
c
               if ( R - rh_o_2 .lt. 0.201d0 ) then
c
                  vlam(ir,1) = 9.999999d+30
                  vlam(ir,2) = 9.999999d+30
                  vlam(ir,3) = 9.999999d+30
c
               else
c
                  if ( isurf .eq. 1 ) then
                     call vH2opt(rh,E,0)
                     Vhh = e(1)
                  else if ( isurf .eq. 2 ) then
                     call vH2opt95(rh,E,0)
                     Vhh = e(1)
                  else if ( isurf .eq. 3 ) then
                     ri(1) = rh
                     ri(2) = 10000.d0
                     ri(3) = 10000.d0
                     call pot
                     Vhh = pote1
                  else if ( isurf .eq. 4 ) then
                     call hhpot(rh,e)
                     Vhh = e(1)
                  else if ( isurf .eq. 5 ) then
                     call h3aguado7(rh,10000.d0,10000.d0,Vhh,dVtot)
                  else if ( isurf .eq. 6 ) then
                     call h3aguado11(rh,10000.d0,10000.d0,Vhh,dVtot)
                  endif
c
                  ntheta = nstart
                  v(1) = 1.d30
                  v(2) = 1.d30
                  v(3) = 1.d30
                  dv = 1.d30
                  dvp = 1.d30
c-debug[
c-debug;               if ( jr .eq. 16 ) write(6,
c-debug;     $              '(/" isurf",i2," r,R,pi",3f15.10/)')
c-debug;     $              isurf, rh, R, pi
c-debug]
c
                  do while ( ntheta .le. nmax .and.
     $                 ( dv .gt. dvtry .or. dvp .gt. dvptry ) )
c
                     dvp = dv
                     vp(1) = v(1)
                     vp(2) = v(2)
                     vp(3) = v(3)
c
                     n2d = nmax / ntheta
                     if ( ntheta .eq. nstart ) then
                        n2b = 0
                        n2e = nmax
                        n2duse = n2d
                     else
                        n2b = n2d
                        n2e = nmax - n2d
                        n2duse = n2d * 2
                     endif
c
                     do i = n2b, n2e, n2duse
c
                        ctheta = i * thninv
                        cthsq = ctheta**2
                        stheta = sqrt( 1.d0 - cthsq )
                        cthetaR = ctheta * R
                        sthetaR = stheta * R
                        ri(1) = rh
                        ri(2) = sqrt( ( cthetaR - rh_o_2 )**2
     $                       + sthetaR**2 )
                        ri(3) = sqrt( ( cthetaR + rh_o_2 )**2
     $                       + sthetaR**2 )
c						   ! pn = (2n+1) * Pn(cosTheta)
                        p2 = 7.5d0 * cthsq - 2.5d0
                        p4 = ( 39.375d0 * cthsq - 33.75d0 )
     $                       * cthsq + 3.375d0
c
                        if ( isurf .eq. 1 ) then
                           call h3tot(ri,Vtot,dVtot,0)
                        else if ( isurf .eq. 2 ) then
                           call bkmp2(ri,Vtot,dVtot,0)
                        else if ( isurf .eq. 3 ) then
                           call pot
                           Vtot = pote1
                        else if ( isurf .eq. 4 ) then
                           call vlsth(ri,Vtot,dvdr1,dvdr2,dvdr3,0,0,0)
                        else if ( isurf .eq. 5 ) then
                           call h3aguado7(ri(1),ri(2),ri(3),Vtot,dVtot)
                        else if ( isurf .eq. 6 ) then
                           call h3aguado11(ri(1),ri(2),ri(3),Vtot,
     $                                     dVtot)
                        endif
                        Vint = Vtot - Vhh
c
                        vi(i,1) = Vint
                        vi(i,2) = Vint * p2
                        vi(i,3) = Vint * p4
c
c-debug[
c-debug;                     if ( jr .eq. 16 .and. ntheta .eq. 8 ) write(6,
c-debug;     $                    '(i7," cos,TH",2f10.6," V,p2,p4",1p,3e24.16)'
c-debug;     $                    ) i, ctheta, 180. * acos(ctheta) / pi,
c-debug;     $                    Vint, p2, p4
c-debug]
                     enddo
c
                     n2duse = n2d * 2
                     do k = 1, 3
                        vv_4 = 0.d0
                        do i2 = n2d, nmax - n2d, n2duse
                           vv_4 = vv_4 + vi(i2,k)
                        enddo
                        vv_2 = 0.d0
                        do i2 = n2duse, nmax - n2duse, n2duse
                           vv_2 = vv_2 + vi(i2,k)
                        enddo
                        v(k) = ( vi(0,k) + vi(nmax,k) + 4.d0 * vv_4
     $                       + 2.d0 * vv_2 ) / ( 3 * ntheta )
                     enddo
c
                     dv = max( abs(v(1)-vp(1))
     $                    / max( abs(v(1)) , abs(vp(1)) , 1.d-6 ) ,
     $                    abs(v(2)-vp(2))
     $                    / max( abs(v(2)) , abs(vp(2)) , 1.d-8 ) ,
     $                    abs(v(3)-vp(3))
     $                    / max( abs(v(3)) , abs(vp(3)) , 1.d-10 ) )
c-debug[
c-debug;                  if ( jr .eq. 16 .and. ntheta .eq. 4 ) write(6,
c-debug;     $                 '(/" ntheta",i7," Vlam,p,d",1p,8e22.14)')
c-debug;     $                 ntheta, v, vp, dv, dvp
c-debug;                  read(5,'(a1)') conf
c-debug]
c
                     ntheta = 2 * ntheta
c
                  enddo
c
                  ntheta = ntheta / 2
c
                  nthmax = max( nthmax , ntheta )
                  dm = max( dm , dv )
                  vlam(ir,1) = v(1)
                  vlam(ir,2) = v(2)
                  vlam(ir,3) = v(3)
c-debug[
c-debug;               if ( jr .eq. 16 ) write(6,
c-debug;     $              '(i2," r,R",2f7.3," n",i7," V,p,d",1p,8e18.10)')
c-debug;     $              isurf, rh, R, ntheta, v, vp, dv, dvp
c-debug;               read(5,'(a1)') conf
c-debug]
c
               endif
c
            enddo
c
            do lam = 0, 4, 2
               k = max( 1 , min( 3 , lam ) )
               write(k,20) R, (vlam(ir,k),ir=1,nri)
 20            format(f5.1,1p,5e14.6)
            enddo
c-debug[
c-debug;            write(6,7347) cend(isurf), R, nthmax, dm,
c-debug;     $           ( lam, ( vlam( ir, max(1,min(3,lam)) ),ir=1,nri),
c-debug;     $           lam = 0, 4, 2 )
c-debug;               if ( jr .eq. 16 .or. jr .eq. 30 .or. jr .eq. 200 )
c-debug;     $              write(6,7347) cend(isurf), R, nthmax, dm, lam,
c-debug;     $              (vlam(ir),ir=1,nri)
c-debug; 7347       format(' ',a5,' R',f5.1,' Nthm',i7,' Dmax',1p,e10.3,
c-debug;     $           ' V',i1,'(ri)',5e14.6/
c-debug;     $           40x,' V',i1,'(ri)',5e14.6/
c-debug;     $           40x,' V',i1,'(ri)',5e14.6)
c-debug]
c
            if ( jr .le. 164 ) then
               jr = jr + 2
            else if ( jr .ge. 170 ) then
               jr = jr + 5
            else
               jr = 170
            endif
c
         enddo
c
         close(1)
         close(2)
         close(3)
c
      enddo
c
      stop
      end
