c
c Compile this program, e.g.:
c
c   f77 -o v0ri_h4bmkp v0ri_h4bmkp.f h4bmkp.f
c
c Execute this program, e.g.:
c
c  ./v0ri_h4bmkp
c
c When it asks, give r1 and r2 in bohrs; note that these distances must not
c  exceed 6 bohrs and their sum r1 + r2 must not exceed 9 bohrs.
c When it asks, give the output filename (blank --> default name); optionally,
c  dVacc (accuracy criterion) and/or Nmax (maximum number of integration points
c  in each dimension) can be specified (on the SAME line, BEFORE the filename).
c When it asks, confirm whether the inputs were correct.
c
c Takes about 5 minutes on [Intel Pentium III, 450 MHz, under Linux]
c  for r1 = r2 = 1.449 bohrs with default dVacc = 1.E-5, Nmax = 128
c  (takes about half an hour with extreme dVacc = 1.E-6, Nmax = 128)
c
c Output file contains the SK and BMKP spherical averages for r1 = r2 = 1.449
c  bohrs, and the BMKP spherical average for the user-supplied r1,r2 values,
c  plus some convergence and anisotropy info.
c
c NOTE: SK surface is defined only for r1 = r2 = 1.449 bohrs; BMKP surface is
c       defined at all r1,r2 but is less accurate for R > 4.5 bohrs or so.
c
      program v0ri_h4bkmp
c     -------------------
c
      implicit double precision (a-h,o-z)
c
c NOTE: nmax (and nstart) must be powers of 2, or program will not work.
c
c-test;      parameter ( dvtry = 1.d-3, nmax = 8, nstart = 2 )
c
c-extremeAcc;      parameter ( dvtry = 1.d-6, nmax = 128, nstart = 2 )
c
      parameter ( dvtry = 1.d-5, nmax = 128, nstart = 2 )
c
      parameter ( c5000= .00000d+0, c6000=-.12137d+2,
     $     c8000=-.22135d+3, c10000=-.50387d+4,
     $     dn000= .46645154723253385d-5, d1000=-.36375d0 )
      parameter ( c5202= .00000d+0, c6202=-.25443d+0,
     $     c8202=-.12176d+2, c10202=-.30724d+3,
     $     dn202= .13041676739095653d-6, d1202=-1.00015d-2 )
      parameter ( c5220= .00000d+0, c6220=-.12569d-1,
     $     c8220=-.23694d+0, c10220=-.12763d+2,
     $     dn220= .51211515143105659d-8, d1220=-.4437d-2 )
      parameter ( c5222= .00000d+0, c6222=-.67187d-2,
     $     c8222= .31517d+0, c10222= .14913d+2,
     $     dn222= .47665894601774043d-9, d1222=2.78035d-3 )
      parameter ( c5224= .12867d+0, c6224=-.40314d-1,
     $     c8224=-.92886d+0, c10224=-.41671d+2,
     $     dn224=-.34612986743414548d-6, d1224=-.8333d-2 )
c
      common/schaei/rschae(91),
     $     v000(91),v202(91),v220(91),v222(91),v224(91),
     $     p000(91),p202(91),p220(91),p222(91),p224(91),
     $     rt5,a220,a222,a224,b220,b222,b224,c220,c222,c224,
     $     ischae
c
      character*255 cfile
c
      character*11 cout(6), ctmp
c
      dimension v1(0:nmax), v1pp(0:nmax), v2(0:nmax), v2p(0:nmax),
     $     v3(0:nmax), cc(4,3), dVdcc(4,3), vout(6)
c
      dimension v000bmkp_1449(91)
c
      data ischae/0/,p000/91*0.d0/,p202/91*0.d0/,p220/91*0.d0/,
     $     p222/91*0.d0/,p224/91*0.d0/
      data rschae/ 2.0,
     $      2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0,
     $      3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0,
     $      4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0,
     $      5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0,
     $      6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0,
     $      7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0,
     $      8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0,
     $      9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9,10.0,
     $     10.1,10.2,10.3,10.4,10.5,10.6,10.7,10.8,10.9,11.0/
      data v000/ .22220d+0,
     $      .18839d+0, .15971d+0, .13537d+0, .11473d+0, .97220d-1,
     $      .82372d-1, .69783d-1, .59110d-1, .50064d-1, .42390d-1,
     $      .35826d-1, .30188d-1, .25378d-1, .21294d-1, .17836d-1,
     $      .14914d-1, .12451d-1, .10378d-1, .86264d-2, .71186d-2,
     $      .58213d-2, .47402d-2, .38528d-2, .31217d-2, .25193d-2,
     $      .20234d-2, .16155d-2, .12811d-2, .10043d-2, .77178d-3,
     $      .57972d-3, .42567d-3, .30276d-3, .20450d-3, .12547d-3,
     $      .63058d-4, .15001d-4,-.22505d-4,-.51408d-4,-.73355d-4,
     $     -.89190d-4,-.99894d-4,-.10657d-3,-.11021d-3,-.11165d-3,
     $     -.11127d-3,-.10943d-3,-.10652d-3,-.10288d-3,-.98782d-4,
     $     -.94356d-4,-.89732d-4,-.85027d-4,-.80332d-4,-.75708d-4,
     $     -.71208d-4,-.66865d-4,-.62703d-4,-.58735d-4,-.54966d-4,
     $     -.51404d-4,-.48052d-4,-.44904d-4,-.41954d-4,-.39194d-4,
     $     -.36617d-4,-.34212d-4,-.31969d-4,-.29881d-4,-.27938d-4,
     $     -.26130d-4,-.24447d-4,-.22880d-4,-.21422d-4,-.20066d-4,
     $     -.18805d-4,-.17632d-4,-.16540d-4,-.15524d-4,-.14577d-4,
     $     -.13694d-4,-.12871d-4,-.12105d-4,-.11390d-4,-.10730d-4,
     $     -.10133d-4,-.95754d-5,-.90481d-5,-.85497d-5,-.80779d-5/
      data v202/ .63286d-2,
     $      .53968d-2, .46017d-2, .39232d-2, .33443d-2, .28505d-2,
     $      .24292d-2, .20700d-2, .17636d-2, .15024d-2, .12796d-2,
     $      .10881d-2, .92246d-3, .77993d-3, .65769d-3, .55321d-3,
     $      .46414d-3, .38842d-3, .32415d-3, .26988d-3, .22490d-3,
     $      .18785d-3, .15677d-3, .13042d-3, .10818d-3, .89435d-4,
     $      .73654d-4, .60392d-4, .49303d-4, .39886d-4, .31687d-4,
     $      .24729d-4, .19063d-4, .14471d-4, .10727d-4, .76356d-5,
     $      .51358d-5, .31712d-5, .17445d-5, .55309d-6,-.45526d-6,
     $     -.12379d-5,-.18068d-5,-.22055d-5,-.24777d-5,-.26584d-5,
     $     -.27595d-5,-.27946d-5,-.27803d-5,-.27297d-5,-.26541d-5,
     $     -.25592d-5,-.24508d-5,-.23339d-5,-.22124d-5,-.20893d-5,
     $     -.19669d-5,-.18468d-5,-.17304d-5,-.16180d-5,-.15093d-5,
     $     -.14059d-5,-.13088d-5,-.12178d-5,-.11327d-5,-.10534d-5,
     $     -.97946d-6,-.91069d-6,-.84680d-6,-.78746d-6,-.73223d-6,
     $     -.68088d-6,-.63332d-6,-.58931d-6,-.54860d-6,-.51096d-6,
     $     -.47617d-6,-.44400d-6,-.41426d-6,-.38674d-6,-.36103d-6,
     $     -.33698d-6,-.31478d-6,-.29433d-6,-.27256d-6,-.25907d-6,
     $     -.24804d-6,-.23866d-6,-.22963d-6,-.22104d-6,-.21227d-6/
      data v220/ .18104d-2,
     $      .14107d-2, .10990d-2, .85606d-3, .66666d-3, .51906d-3,
     $      .40406d-3, .31447d-3, .24470d-3, .19038d-3, .14803d-3,
     $      .11467d-3, .87983d-4, .66725d-4, .49854d-4, .36533d-4,
     $      .26083d-4, .17951d-4, .11607d-4, .69747d-5, .46418d-5,
     $      .40932d-5, .38530d-5, .34351d-5, .30137d-5, .26104d-5,
     $      .22371d-5, .18997d-5, .16019d-5, .13335d-5, .10795d-5,
     $      .84995d-6, .65872d-6, .49861d-6, .35861d-6, .22336d-6,
     $      .10544d-6, .16344d-7,-.51097d-7,-.10950d-6,-.17257d-6,
     $     -.22351d-6,-.25597d-6,-.27366d-6,-.28713d-6,-.31459d-6,
     $     -.33452d-6,-.33810d-6,-.33139d-6,-.32450d-6,-.33855d-6,
     $     -.34981d-6,-.34620d-6,-.33399d-6,-.31609d-6,-.29464d-6,
     $     -.27135d-6,-.24738d-6,-.22283d-6,-.20309d-6,-.19821d-6,
     $     -.19698d-6,-.19042d-6,-.18074d-6,-.16911d-6,-.15642d-6,
     $     -.14329d-6,-.13019d-6,-.11704d-6,-.10534d-6,-.10097d-6,
     $     -.99782d-7,-.96220d-7,-.91266d-7,-.85492d-7,-.79308d-7,
     $     -.73001d-7,-.65024d-7,-.56352d-7,-.48175d-7,-.41179d-7,
     $     -.35195d-7,-.30077d-7,-.25700d-7,-.21971d-7,-.19008d-7,
     $     -.16391d-7,-.13943d-7,-.11891d-7,-.10121d-7,-.86923d-8/
      data v222/-.11466d-2,
     $     -.89589d-3,-.69983d-3,-.54656d-3,-.42678d-3,-.33318d-3,
     $     -.26005d-3,-.20294d-3,-.15833d-3,-.12352d-3,-.96311d-4,
     $     -.75152d-4,-.58604d-4,-.45614d-4,-.35389d-4,-.27333d-4,
     $     -.20989d-4,-.16003d-4,-.12053d-4,-.90829d-5,-.74079d-5,
     $     -.67332d-5,-.62384d-5,-.56592d-5,-.50959d-5,-.45624d-5,
     $     -.40664d-5,-.36114d-5,-.32003d-5,-.28217d-5,-.24618d-5,
     $     -.21341d-5,-.18518d-5,-.16083d-5,-.13953d-5,-.12031d-5,
     $     -.10346d-5,-.89267d-6,-.81799d-6,-.72885d-6,-.62326d-6,
     $     -.52750d-6,-.44948d-6,-.38627d-6,-.33160d-6,-.27816d-6,
     $     -.23250d-6,-.19695d-6,-.16885d-6,-.14467d-6,-.11787d-6,
     $     -.94557d-7,-.77707d-7,-.65258d-7,-.56131d-7,-.49511d-7,
     $     -.44720d-7,-.41270d-7,-.38983d-7,-.36304d-7,-.30477d-7,
     $     -.24345d-7,-.20216d-7,-.17447d-7,-.15679d-7,-.14631d-7,
     $     -.14092d-7,-.13901d-7,-.14002d-7,-.14096d-7,-.13104d-7,
     $     -.11608d-7,-.10507d-7,-.95111d-8,-.86068d-8,-.77865d-8,
     $     -.70427d-8,-.63685d-8,-.58042d-8,-.52504d-8,-.47485d-8,
     $     -.42936d-8,-.38817d-8,-.35087d-8,-.31724d-8,-.28927d-8,
     $     -.26303d-8,-.23704d-8,-.21389d-8,-.19268d-8,-.17473d-8/
      data v224/ .45157d-2,
     $      .37475d-2, .31095d-2, .25798d-2, .21400d-2, .17748d-2,
     $      .14718d-2, .12203d-2, .10116d-2, .83854d-3, .69491d-3,
     $      .57707d-3, .48075d-3, .40170d-3, .33657d-3, .28278d-3,
     $      .23824d-3, .20129d-3, .17049d-3, .14507d-3, .12515d-3,
     $      .10972d-3, .96767d-4, .85413d-4, .75562d-4, .66999d-4,
     $      .59542d-4, .53036d-4, .47352d-4, .42364d-4, .37972d-4,
     $      .34110d-4, .30720d-4, .27736d-4, .25099d-4, .22761d-4,
     $      .20703d-4, .18905d-4, .18320d-4, .17293d-4, .15895d-4,
     $      .14542d-4, .13354d-4, .12311d-4, .11368d-4, .10511d-4,
     $      .97315d-5, .90231d-5, .83777d-5, .77877d-5, .72454d-5,
     $      .67485d-5, .62936d-5, .58761d-5, .54923d-5, .51390d-5,
     $      .48130d-5, .45120d-5, .42335d-5, .39755d-5, .37364d-5,
     $      .35148d-5, .33091d-5, .31178d-5, .29398d-5, .27738d-5,
     $      .26190d-5, .24743d-5, .23392d-5, .22126d-5, .20933d-5,
     $      .19814d-5, .18772d-5, .17797d-5, .16885d-5, .16029d-5,
     $      .15226d-5, .14470d-5, .13759d-5, .13088d-5, .12455d-5,
     $      .11856d-5, .11290d-5, .10757d-5, .10254d-5, .97779d-6,
     $      .93230d-6, .88884d-6, .84737d-6, .80777d-6, .77024d-6/
c
      data v000bmkp_1449 / 2.1129d-1,
     $      1.8023d-1,  1.5389d-1,  1.3144d-1,  1.1226d-1,  9.5828d-2,
     $      8.1743d-2,  6.9653d-2,  5.9268d-2,  5.0340d-2,  4.2665d-2,
     $      3.6068d-2,  3.0406d-2,  2.5554d-2,  2.1408d-2,  1.7876d-2,
     $      1.4879d-2,  1.2347d-2,  1.0216d-2,  8.4323d-3,  6.9444d-3,
     $      5.7074d-3,  4.6815d-3,  3.8315d-3,  3.1272d-3,  2.5430d-3,
     $      2.0575d-3,  1.6535d-3,  1.3170d-3,  1.0368d-3,  8.0401d-4,
     $      6.1142d-4,  4.5310d-4,  3.2406d-4,  2.1998d-4,  1.3705d-4,
     $      7.1863d-5,  2.1390d-5, -1.7054d-5, -4.5805d-5, -6.6859d-5,
     $     -8.1878d-5, -9.2219d-5, -9.8963d-5, -1.0296d-4, -1.0485d-4,
     $     -1.0514d-4, -1.0420d-4, -1.0231d-4, -9.9704d-5, -9.6539d-5,
     $     -9.2958d-5, -8.9072d-5, -8.4973d-5, -8.0740d-5, -7.6436d-5,
     $     -7.2117d-5, -6.7829d-5, -6.3608d-5, -5.9486d-5, -5.5487d-5,
     $     -5.1630d-5, -4.7930d-5, -4.4398d-5, -4.1039d-5, -3.7858d-5,
     $     -3.4856d-5, -3.2033d-5, -2.9385d-5, -2.6908d-5, -2.4599d-5,
     $     -2.2450E-5, -2.0455E-5, -1.8606E-5, -1.6897E-5, -1.5321E-5,
     $     -1.3869E-5, -1.2534E-5, -1.1309E-5, -1.0187E-5, -9.1607E-6,
     $     -8.2243E-6, -7.3712E-6, -6.5953E-6, -5.8908E-6, -5.2519E-6,
     $     -4.6731E-6, -4.1492E-6, -3.6755E-6, -3.2476E-6, -2.8615E-6 /
c
      pi = acos(-1.d0)
c
      n_max = nmax
      d_vtry = dvtry
c
      fri = 0.d0
      rt5=dsqrt(5.d0)
      c220=rt5
      b220=-3.d0*rt5
      a220=.75d0*rt5
      c222=-5.d0/dsqrt(0.7d0)
      b222=-1.5d0*c222
      a222=-.75d0*c222
      c224=45.d0/dsqrt(17.5d0)
      b224=2.d0*c224
      a224=.125d0*c224
c
c  Formula below gives V_SK at given R, theta1, theta2, phi at r1=r2=1.449 bohr
c
c-SK;      y1=1.5d0*costh1**2-.5d0
c-SK;      y2=1.5d0*costh2**2-.5d0
c-SK;      vsk=vr000 + rt5*vr202*(y1+y2)
c-SK;     $     + ( (a220*vr220+a222*vr222+a224*vr224)
c-SK;     $     *(sinth1*sinth2)**2*cos2ph
c-SK;     $     + (b220*vr220+b222*vr222+b224*vr224)
c-SK;     $     *costh1*costh2*sinth1*sinth2*cosphi
c-SK;     $     + (c220*vr220+c222*vr222+c224*vr224)*y1*y2 )
c
      h202 = rt5
      h220 = max( abs(a220) , abs(b220) , abs(c220) )
      h222 = max( abs(a222) , abs(b222) , abs(c222) )
      h224 = max( abs(a224) , abs(b224) , abs(c224) )
c
c-SK;      r1 = 1.449d0
c-SK;      r2 = 1.449d0
c
      write(6,2)
 2    format(/' Give  r1 r2: ',$)
      read(5,*) r1, r2
      if ( min(r1,r2) .lt. 0.001d0 .or. r1 + r2 .gt. 9.000001d0 .or.
     $     max(r1,r2) .gt. 6.000001d0 ) stop
     $     ' STOP -- Error: bad input r1,r2 '
c
      write(6,4)
 4    format(/' Give  [ dVacc(opt) ] [ Nmax(opt) ] output_file_name:')
      read(5,'(a155)') cfile
      if ( cfile(255:255) .ne. ' ' ) stop
     $     ' STOP -- Error: output filename too long. '
c
      last = 255
      do while ( last .gt. 1 .and. cfile(last:last) .eq. ' ' )
         last = last - 1
      enddo
      k = 1
      do while ( k .lt. last .and. cfile(k:k) .eq. ' ' )
         k = k + 1
      enddo
c
      if ( cfile(k:k) .ge. '0' .and. cfile(k:k) .le. '9' ) then
c
         k0 = k
         n = 0
         do while ( n .lt. 100000000 .and.
     $        cfile(k:k) .ge. '0' .and. cfile(k:k) .le. '9' )
            read(cfile(k:k),'(i1)') i
            n = 10 * n + i
            k = k + 1
         enddo
c
         if ( cfile(k:k) .eq. ' ' ) then
c
            n_max = min( nmax , max( n , 2 * nstart ) )
            k0 = k
c
         else if ( cfile(k:k) .eq. '.' ) then
c
            v = n
            den = 1.d0
            k = k + 1
            do while ( cfile(k:k) .ge. '0' .and. cfile(k:k) .le. '9' )
               den = den / 10.d0
               read(cfile(k:k),'(i1)') i
               v = v + i / den
               k = k + 1
            enddo
            if ( cfile(k:k) .eq. ' ' ) then
               d_vtry = max( 0.d0 , min( 100.d0 , v ) )
               k0 = k
            else if ( cfile(k:k) .eq. 'e' .or. cfile(k:k) .eq. 'E' .or.
     $              cfile(k:k) .eq. 'd' .or. cfile(k:k) .eq. 'D' ) then
               k = k + 1
               if ( cfile(k:k) .eq. '-' ) then
                  isgn = -1
                  k = k + 1
               else
                  isgn = 1
                  if ( cfile(k:k) .eq. '+' ) k = k + 1
               endif
               n = 0
               do while ( n .lt. 1000 .and.
     $              cfile(k:k) .ge. '0' .and. cfile(k:k) .le. '9' )
                  read(cfile(k:k),'(i1)') i
                  n = 10 * n + i
                  k = k + 1
               enddo
               if ( cfile(k:k) .eq. ' ' ) then
                  d_vtry = max( 0.d0 , min( 100.d0 ,
     $                 v * 10.d0**( isgn * n ) ) )
                  k0 = k
               endif
            endif
c
            if ( k0 .eq. k ) then
               do while ( k .lt. last .and. cfile(k:k) .eq. ' ' )
                  k = k + 1
               enddo
               k0 = k
               if ( cfile(k:k) .ge. '0' .and.
     $              cfile(k:k) .le. '9' ) then
                  n = 0
                  do while ( n .lt. 100000000 .and.
     $                 cfile(k:k) .ge. '0' .and. cfile(k:k) .le. '9' )
                     read(cfile(k:k),'(i1)') i
                     n = 10 * n + i
                     k = k + 1
                  enddo
                  if ( cfile(k:k) .eq. ' ' ) then
                     n_max = min( nmax , max( n , 2 * nstart ) )
                     k0 = k
                  endif
               endif
            endif
c
         endif
c
         k = k0
         do while ( k .lt. last .and. cfile(k:k) .eq. ' ' )
            k = k + 1
         enddo
c
      endif
c
      if ( k .gt. last .or. cfile .eq. ' ' ) then
         write(cfile,6) r1, r2
 6       format('v0ri_h4_',f8.6,'_',f8.6,'_BMKP.txt')
         last = 255
      else if ( k .gt. 1 ) then
         j = 1
         do i = k, last
            cfile(j:j) = cfile(i:i)
            j = j + 1
         enddo
         cfile(j:) = ' '
      endif
c
      do while ( last .gt. 1 .and. cfile(last:last) .eq. ' ' )
         last = last - 1
      enddo
c
      do k = 1, last
         if ( cfile(k:k) .eq. ' ' ) cfile(k:k) = '_'
      enddo
c
      n = nmax
      do while ( n / 2 .ge. max( n_max , 2 * nstart ) )
         n = n / 2
      enddo
      n_max = n
      d_vtry = max( 0.d0 , min( 100.d0 , d_vtry ) )
c
      write(6,8) r1, r2, d_vtry, n_max, cfile(:last)
 8    format(/' r1=',f13.6,'  r2=',f13.6,'     [ dVacc=',1p,e9.2,
     $     '  Nmax=',i5,' ]'/' outfile= ',a/' Confirm? (y/n): ',$)
      read(5,'(a1)') ctmp(1:1)
c
      if ( ctmp(1:1) .ne. 'y' .and. ctmp(1:1) .ne. 'Y' ) stop
     $     ' STOP -- Run cancelled. '
c
      vinvnmax = 1.d0 / n_max
      r1o2 = r1 / 2.d0
      r2o2 = r2 / 2.d0
c
      call vH2opt95_h4bmkp( r1, vout, 0 )
      call vH2opt95_h4bmkp( r2, vout(2), 0 )
      v2hh = vout(1) + vout(2)
c
      do i = 1, 4
         cc(i,1) = 0.d0
         cc(i,2) = 0.d0
         cc(i,3) = 0.d0
      enddo
c
      open(1,file=cfile,form='FORMATTED',status='UNKNOWN')
c
c-debug[
c-debug;      write(6,10) d_vtry, n_max, r1, r2
c-debug]
      write(1,10) d_vtry, n_max, r1, r2
 10   format('# Spherical average V_000 of the H4',
     $     ' surface for H2-size r1 = r2 = 1.449 bohrs'/
     $     '# (SK and BMKP surfaces) and for USER-SUPPLIED',
     $     ' {r1,r2} (for BMKP surface only);'/
     $     '# energies in hartrees, distances in bohrs.'/
     $     '#   SK = J. Schaefer and W.E. Kohler,',
     $     ' Z. Phys. D 13, 217 (1989).'/
     $     '#   BMKP = A.I. Boothroyd, P.G. Martin,',
     $     ' W.J. Keogh, and M.R. Peterson,'/
     $     '#           J. Chem. Phys., 116, 666-689 (2002).'/
     $     '# Note that the BMKP surface should be',
     $     ' more accurate for R < 5 bohrs, roughly.'/
     $     '#'/
     $     '# For the SK surface, the absolute value',
     $     ' of the largest anisotropic term is'/
     $     '#   also given (max|Vani| as well as V_000).'/
     $     '# For the BMKP surface, V_000 is given both',
     $     ' for r1 = r2 = 1.449 and for the'/
     $     '#   USER-SUPPLIED {r1,r2}; these were integrated',
     $     ' numerically, and are only'/
     $     '#   calculated out to R = 11.0 bohrs.  V[n/8] is',
     $     ' the numerical integral for a'/
     $     '#   case where each of the 3 dimensions of the',
     $     ' integral used only half as many'/
     $     '#   points as the final integral did; |dV[n/8]|',
     $     ' is the absolute FRACTIONAL'/
     $     '#   difference between V and V[n/8].  |dVhilo|',
     $     ' gives the absolute FRACTIONAL'/
     $     '#   difference between the maximum and minimum',
     $     ' values of the BMKP H4 surface'/
     $     '#   over all points used in the integral',
     $     ' (a measure of the BMKP anisotropy).'/
     $     '#'/
     $     '#     ==========........... ========== **********',
     $     '   dVacc=',1p,e9.2,'  Nmax=',i5/
     $     '#   r1= 1.449      1.449      1.449 ',0p,f13.6/
     $     '#   r2= 1.449      1.449      1.449 ',0p,f13.6,
     $     '           ....fractional....'/
     $     '#  R   V_000(SK) max|aVan| V_000(BMKP) V_000(BMKP',
     $     ')   V[n/8]  |dV[n/8]| |dVhilo|'/
     $     '#==== ========== ---------- ========== **********',
     $     ' ---------- -------- ---------')
c
c Spherical average V_000 of the H4 surface for H2-size r1 = r2 = 1.449 bohrs
c (SK and BMKP surfaces) and for USER-SUPPLIED {r1,r2} (for BMKP surface only);
c energies in hartrees, distances in bohrs.
c   SK = J. Schaefer and W.E. Kohler, Z. Phys. D 13, 217 (1989).
c   BMKP = A.I. Boothroyd, P.G. Martin, W.J. Keogh, and M.R. Peterson,
c           J. Chem. Phys., 116, 666-689 (2002).
c Note that the BMKP surface should be more accurate for R < 5 bohrs, roughly.
c
c For the SK surface, the absolute value of the largest anisotropic term is
c   also given (max|aVan| as well as V_000).
c For the BMKP surface, V_000 is given both for r1 = r2 = 1.449 and for the
c   USER-SUPPLIED {r1,r2}; these were integrated numerically, and are only
c   calculated out to R = 11.0 bohrs.  V[n/8] is the numerical integral for a
c   case where each of the 3 dimensions of the integral used only half as many
c   points as the final integral did; |dV[n/8]| is the absolute FRACTIONAL
c   difference between V and V[n/8].  |dVhilo| gives the absolute FRACTIONAL
c   difference between the maximum and minimum values of the BMKP H4 surface
c   over all points used in the integral (a measure of the BMKP anisotropy).
c
c     ==========........... ========== **********   dVacc= 1.00E-05  Nmax=  128
c   r1= 1.449      1.449      1.449      1.200000
c   r2= 1.449      1.449      1.449      1.600000           ....fractional....
c  R   V_000(SK) max|aVan| V_000(BMKP) V_000(BMKP)   V[n/8]  |dV[n/8]| |dVhilo|
c==== ========== ---------- ========== ********** ---------- -------- ---------
c 9.0 -1.8839E+0 -1.8839E+0 -1.8839E+0 -1.8839E+0 -1.8839E+0 1.31E-07 1.644E-01
c10.1 -1.8839E+0 -1.8839E+0  0.0000E+0 -1.8839E+0
c
c
c  Calculate integral, for R from 2 bohrs to 11 bohrs
c
      do jr = 1, 91
c
         R = rschae(jr)
c
         if ( R - r1o2 - r2o2 .lt. 0.001d0 ) then
c
            v = 9.999d+30
            vppp = 0.d0
            dvm = 1.d0
            dvno8 = 1.d0
c
            n1 = 0
            n2hi = 0
            n3hi = 0
c
         else
c
            v_max = -1.d30
            v_min = 1.d30
            dv = 1.d30
            v = 1.d30
c
            n1 = nstart
            n2hi = 0
            n3hi = 1
c
            do while ( n1 .le. n_max .and. dv .gt. d_vtry )
c
               vp = v
c
               n1d = n_max / n1
               if ( n1 .eq. nstart ) then
                  n1b = 0
                  n1e = n_max
                  n1duse = n1d
               else
                  n1b = n1d
                  n1e = n_max - n1d
                  n1duse = n1d * 2
               endif
c
               do i1 = n1b, n1e, n1duse
c
                  cth1 = min( 1.d0 , i1 * vinvnmax )
                  sth1 = sqrt( 1.d0 - cth1**2 )
                  cc(2,2) = cth1 * r1o2 - R
                  cc(2,3) = sth1 * r1o2
                  cc(1,2) = - cth1 * r1o2 - R
                  cc(1,3) = - cc(2,3)
c
                  dv1 = 1.d30
                  v1(i1) = 1.d30
c
                  n2 = nstart
c
                  do while ( n2 .le. n_max .and. dv1 .gt. d_vtry )
c
                     v1p = v1(i1)
c
                     n2d = n_max / n2
                     if ( n2 .eq. nstart ) then
                        n2b = 0
                        n2e = n_max
                        n2duse = n2d
                     else
                        n2b = n2d
                        n2e = n_max - n2d
                        n2duse = n2d * 2
                     endif
c
                     do i2 = n2b, n2e, n2duse
c
                        cth2 = min( 1.d0 , i2 * vinvnmax )
                        sth2 = sqrt( 1.d0 - cth2**2 )
                        cc(4,2) = r2o2 * cth2
                        cc(3,2) = - cc(4,2)
                        rho2 = r2o2 * sth2
c
                        dv2 = 1.d30
                        v2(i2) = 1.d30
c
                        if ( i1 .eq. n_max .or. i2 .eq. n_max ) then
                           nmax_p = 0
                        else
                           nmax_p = n_max
                        endif
c
                        n3 = min( nstart , nmax_p )
c
                        do while ( n3 .le. nmax_p .and.
     $                       dv2 .gt. d_vtry )
c
                           v2p(i2) = v2(i2)
c
                           if ( nmax_p .eq. 0 ) then
                              n3d = 1
                              n3b = 0
                              n3e = 0
                              n3duse = 1
                           else if ( n3 .eq. nstart ) then
                              n3d = nmax_p / n3
                              n3b = 0
                              n3e = nmax_p
                              n3duse = n3d
                           else
                              n3d = nmax_p / n3
                              n3b = n3d
                              n3e = nmax_p - n3d
                              n3duse = n3d * 2
                           endif
c
                           phi = 0.d0
c
                           do i3 = n3b, n3e, n3duse
c
                              if ( nmax_p .gt. 0 )
     $                             phi = ( pi * i3 ) / nmax_p
                              cc(3,3) = rho2 * cos(phi)
                              cc(3,1) = rho2 * sin(phi)
                              cc(4,3) = - cc(3,3)
                              cc(4,1) = - cc(3,1)
c
                              call h4bmkp_cc( cc, vtot, dVdcc, 0 )
                              v3(i3) = vtot - v2hh
                              v_max = max( v_max , v3(i3) )
                              v_min = min( v_min , v3(i3) )
c
                           enddo
c
                           if ( nmax_p .eq. 0 ) then
                              v2(i2) = v3(0)
                              v2p(i2) = v3(0)
                              dv2 = 0.d0
                           else
                              n3duse = n3d * 2
                              vv_4 = 0.d0
                              do i3 = n3d, nmax_p - n3d, n3duse
                                 vv_4 = vv_4 + v3(i3)
                              enddo
                              vv_2 = 0.d0
                              do i3 = n3duse, nmax_p - n3duse, n3duse
                                 vv_2 = vv_2 + v3(i3)
                              enddo
                              v2(i2) = ( v3(0) + v3(nmax_p) + 4.d0
     $                             * vv_4 + 2.d0 * vv_2 ) / ( 3 * n3 )
                              dv2 = abs( v2(i2) - v2p(i2) )
     $                             / max( 1.d-5 , abs( v2(i2) ) )
                           endif
c-debug[
c-debug;                        write(6,'("  n3=",i3," v2(",i3,")",1p,e14.6,
c-debug;     $                       " <--",9e14.6)')
c-debug;     $                       n3, i2, v2(i2), (v3(i3),i3=0,nmax_p,n3d)
c-debug]
c
                           n3 = max( 1 , 2 * n3 )
c
                        enddo
c
                        n3hi = max( n3hi , n3 / 2 )
c
                     enddo
c
                     n2duse = n2d * 2
                     vv_4 = 0.d0
                     do i2 = n2d, n_max - n2d, n2duse
                        vv_4 = vv_4 + v2(i2)
                     enddo
                     vv_2 = 0.d0
                     do i2 = n2duse, n_max - n2duse, n2duse
                        vv_2 = vv_2 + v2(i2)
                     enddo
                     v1(i1) = ( v2(0) + v2(n_max) + 4.d0 * vv_4
     $                    + 2.d0 * vv_2 ) / ( 3 * n2 )
                     dv1 = abs( v1(i1) - v1p )
     $                    / max( 1.d-5 , abs( v1(i1) ) )
c-debug[
c-debug;                  write(6,'(" n2=",i3," v1(",i3,")",1p,e14.6,
c-debug;     $                 " <--",9e14.6)')
c-debug;     $                 n2, i1, v1(i1), (v2(i2),i2=0,n_max,n2d)
c-debug]
c
                     n2 = 2 * n2
c
                  enddo
c
                  n2hi = max( n2hi , n2 / 2 )
c
                  n2 = n2 / 4
                  n2d = n_max / n2
                  n2duse = 2 * n2d
                  vv_4 = 0.d0
                  do i2 = n2d, n_max - n2d, n2duse
                     vv_4 = vv_4 + v2p(i2)
                  enddo
                  vv_2 = 0.d0
                  do i2 = n2duse, n_max - n2duse, n2duse
                     vv_2 = vv_2 + v2p(i2)
                  enddo
                  v1pp(i1) = ( v2p(0) + v2p(n_max) + 4.d0 * vv_4
     $                 + 2.d0 * vv_2 ) / ( 3 * n2 )
c-debug[
c-debug;               write(6,'(" /2=",i3," v1pp( )",1p,e14.6,
c-debug;     $              " <--",9e14.6)')
c-debug;     $              n2, v1pp(i1), (v2p(i2),i2=0,n_max,n2d)
c-debug]
c
               enddo
c
               n1duse = 2 * n1d
               vv_4 = 0.d0
               do i1 = n1d, n_max - n1d, n1duse
                  vv_4 = vv_4 + v1(i1)
               enddo
               vv_2 = 0.d0
               do i1 = n1duse, n_max - n1duse, n1duse
                  vv_2 = vv_2 + v1(i1)
               enddo
               v = ( v1(0) + v1(n_max) + 4.d0 * vv_4 + 2.d0 * vv_2 )
     $              / ( 3 * n1 )
               dv = abs( v - vp ) / max( 1.d-5 , abs(v) )
c-debug[
c-debug;            write(6,'("n1=",i3," v......",1p,e14.6,
c-debug;     $           " <--",9e14.6)')
c-debug;     $           n1, v, (v1(i1),i1=0,n_max,n1d)
c-debug]
c
               n1 = 2 * n1
c
            enddo
c
            n1 = n1 / 4
            n1d = n_max / n1
            n1duse = 2 * n1d
            vv_4 = 0.d0
            do i1 = n1d, n_max - n1d, n1duse
               vv_4 = vv_4 + v1pp(i1)
            enddo
            vv_2 = 0.d0
            do i1 = n1duse, n_max - n1duse, n1duse
               vv_2 = vv_2 + v1pp(i1)
            enddo
            vppp = ( v1pp(0) + v1pp(n_max) + 4.d0 * vv_4
     $           + 2.d0 * vv_2 ) / ( 3 * n1 )
c-debug[
c-debug;         write(6,'("/2=",i3," vppp...",1p,e14.6,
c-debug;     $        " <--",9e14.6)')
c-debug;     $        n1, vppp, (v1pp(i1),i1=0,n_max,n1d)
c-debug]
c
            dvm = abs( v_max - v_min )
     $           / max( 1.d-6 , abs(v_max) , abs(v_min) )
            dvno8 = abs( v - vppp ) / max( 1.d-6 , abs(v) )
c
         endif
c
         vout(1) = v000(jr)
         vout(2) = max( abs( h202 * v202(jr) ) ,
     $        abs( h220 * v220(jr) ) , abs( h222 * v222(jr) ) ,
     $        abs( h224 * v224(jr) ) )
         vout(3) = v000bmkp_1449(jr)
         vout(4) = v
         vout(5) = vppp
c
         do i = 1, 5
            write(ctmp,'(1p,e11.4)') vout(i)
            if ( ctmp(8:8) .ne. '+' .and. ctmp(8:8) .ne. '-' .and.
     $           ctmp(10:10) .eq. '0' ) then
               cout(i) = ' ' // ctmp(1:9) // ctmp(11:11)
            else
               write(cout(i),'(1p,e11.3)') vout(i)
            endif
         enddo
c
         write(1,20) R, ( cout(i), i = 1, 5 ), dvno8, dvm
 20      format(f5.1,5a11,1p,e9.2,e10.3)
c-debug[
c-debug;         write(6,'(" ")')
c-debug;         write(6,'(" n1,23hi",3i4,2(" ",a,1p,e13.5)," dv",3e8.1)')
c-debug;     $        2 * n1, n2hi, n3hi, 'v', v, 'vppp', vppp, dv, dvno8,
c-debug;     $        abs(v/v000bmkp_1449(jr)-1.d0)
c-debug;         write(6,'(" ")')
c-debug]
c-debug[
c-debug;         write(6,20) R, ( cout(i), i = 1, 5 ), dvno8, dvm
c-debug]
c-debug[
c-debug;         write(6,'(" ")')
c-debug]
c
      enddo
c
c  R > 11 bohrs
c
      jr = 112
c
      do while ( jr .le. 200 )
c
         R = 0.1d0 * jr
         r2=r*r
         r6=r2*r2*r2
         vout(1) = ((c10000/r2+c8000)/r2+c6000)/r6
         vout(2) = ((c10202/r2+c8202)/r2+c6202)/r6
         vout(3) = ((c10220/r2+c8220)/r2+c6220)/r6
         vout(4) = ((c10222/r2+c8222)/r2+c6222)/r6
         vout(5) = ((c10224/r2+c8224)/r2+c6224+c5224*r)/r6
         vout(2) = max( abs( h202 * vout(2) ) ,
     $        abs( h220 * vout(3) ) , abs( h222 * vout(4) ) ,
     $        abs( h224 * vout(5) ) )
c
         do i = 1, 2
            write(ctmp,'(1p,e11.4)') vout(i)
            if ( ctmp(10:10) .eq. '0' ) then
               cout(i) = ' ' // ctmp(1:9) // ctmp(11:11)
            else
               write(cout(i),'(1p,e11.3)') vout(i)
            endif
         enddo
c
         write(1,20) R, ( cout(i), i = 1, 2 )
c
         if ( jr .lt. 120 ) then
            jr = jr + 2
         else
            jr = jr + 5
         endif
c
      enddo
c
      close(1)
c
      stop
      end

