c
c  _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
c  _/                                                                _/
c  _/        He-H2 surface evaluation routines (BMP surface)         _/
c  _/        (from the fit called "fit5991.far3.nf12")               _/
c  _/        ---Version of 12 June 2002: Arnold I. Boothroyd         _/
c  _/                                                                _/
c  _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
c
c
c  Any QUESTIONS/PROBLEMS/COMMENTS concerning this programme can be
c  addressed to: boothroy@cita.utoronto.ca or pgmartin@cita.utoronto.ca
c
c
c  NOTE: all real variables are double precision (real*8).
c
c  Distances are in bohrs, energies in hartrees.
c  Note: 1 bohr = 0.52917725 Angstroms, 1 hartree = 627.5096 kcal/mole.
c
c
c  This file contains the subroutines (described in more detail below):
c  --------------------------------------------------------------------
c  heh2bmp_fix_close, heh2bmp_set_err, heh2bmp_ask_set, heh2bmp_cc,
c  heh2bmp_cc3, heh2bmp_ang, heh2bmp, heh2MR, VintMR, heh2modMR,
c  VintmodMR, heh2MR_cc, VintMR_cc, heh2MR_cc3, VintMR_cc3, heh2MR_ang,
c  VintMR_ang, heh2MR_nD, VintMR_nD, Vhh_heh2bmp, geomchk_heh2bmp
c
c
c  The BMP He-H2 surface is always continuous with continuous first
c  derivatives.  The Vhh curve has discontinuous second derivative at
c  r = 0.0451564988896124 (and discontinuous third derivative at
c  r = 0.1193892052031558); in some cases, there may be other places
c  where the second (or third) derivative of the interaction energy is
c  discontinuous, as described below.  The interaction region is well
c  fitted, except for the "rounding off" of the conical intersection
c  of the ground state with the first excited state (at energies above
c  twice the H2-dissociation energy, relative to H2 + He).  The H2 - He
c  van der Waals well is accurately fitted, as is the asymptotic (long
c  range) H2 - He interaction.
c
c  r is the distance between the two H atoms
c  Ra is the distance between the second H atom and the He atom
c  Rb is the distance between the first H atom and the He atom
c  RR is the distance between the He atom and the H2 molecule (center)
c  g is the angle between r and RR
c  a2, y3, z3 are the minimal cartesian coordinates (that were used
c               when specifying the points fitted by the BMP surface)
c
c                        z-axis         y = y3
c                       ^               |
c                       |               |
c                       |               V
c                       |         Ra
c                       |    ______-----He  <-- z = z3
c      ---   z = a2 --> H----         _// 3.
c       |               | 2.        _/ /
c       |               |         _/  /
c       |               |       _/   /
c       |               |__g  _/    /
c       |               |  \_/ RR  / Rb
c       |               | _/      /
c       |               |/       /                     y-axis
c     r |           ----+-------/--------------------->
c       |               |      /
c       |               |     /
c       |               |    /
c       |               |   /
c       |               |  /
c       |               | /
c       |               |/
c      ---  z = -a2 --> H
c                       | 1.
c                       |
c
c######################################################################
c
c  The extrapolation of the BMP He-H2 surface to short distances can be
c  controlled by the user by setting the (integer) close-extrapolation
c  flag  i_fix_close  via:
c
c     call heh2bmp_fix_close( i_fix_close )
c
c  i_fix_close = 0 (DEFAULT): the fitted BMP He-H2 surface rises very
c                             steeply for Ra or Rb < 0.6 bohrs.  If
c                             such a hard repulsive H-He core is not
c                             acceptable, use one of the alternatives
c                             below (i_fix_close = 2 is best).
c
c  i_fix_close = 1 : Switch over to the modMR (modified Muchnick and
c                    Russek) surface (or an exponentially damped 2/R
c                    potential at large r) as Ra or Rb is reduced
c                    from 0.7 to 0.6 bohr.  The third derivative is
c                    discontinuous at Ra or Rb = 0.6, 0.65, 0.7 bohr.
c                 *** NOT RECOMMENDED: there can be a significant
c                    "wiggle" in the surface between 0.6 and 0.7 bohr.
c
c  i_fix_close = 2 : As Ra or Rb drops from 0.65 to 0.55 bohr, switch
c                    over to a distance-dependent factor times the
c                    value of the surface at larger Ra,Rb; the second
c                    derivative is discontinuous at Ra or Rb = 0.55 or
c                    0.65 bohr.  This is the best of the alternatives
c                    to i_fix_close = 0, although the angular variation
c                    of the surface at small r is a bit excessive.
c
c  i_fix_close = 3 : The same as i_fix_close = 2, except that there is
c                    an additional switchover to the modMR surface as
c                    r is reduced from 0.5 to 0.3 bohr.  The second
c                    derivative is discontinuous at Ra or Rb = 0.55 or
c                    0.65 bohr; the third derivative is discontinuous
c                    at r = 0.3, 0.4, and 0.5 bohr.
c                 *** NOT RECOMMENDED: the angular variation of the
c                    surface at small r is reduced, but at the expense
c                    of a significant "wiggle" in the surface between
c                    r = 0.3 and 0.5 bohr.
c
c######################################################################
c
c  The behavior of the program when an error is encountered, and the
c  definition of what the program considers an error, can be controlled
c  by the user by setting the (integer) error level via:
c
c     call heh2bmp_set_err( level )
c
c  level = 3 : The validity of the geometry is checked (i.e., it is an
c              error if the sum of the shortest two distances is less
c              than 1.E-6 bohr shorter than the longest distance).  It
c              is also an error if r < 0.4 bohr or Ra or Rb < 0.6 bohr.
c              On error, a value of V_err = 1.0E+36 is returned for V,
c              Vhh, and Vint, and the derivatives are set to zero.
c
c  level = 2 : Geometry validity is checked; it is also an error if
c              r < 0.001 bohr or Ra or Rb < 0.6 bohr.
c
c  level = 1 (DEFAULT) : Geometry validity is checked; it is also an
c                        error if r or Ra or Rb < 0.001 bohr.
c
c  level = 0 : No error checking is performed (NOT RECOMMENDED: very
c              short distances can cause overflow errors).
c
c  level = -1 to -3 : same as 1 to 3 respectively, except that an error
c                     causes the program to stop.
c
c######################################################################
c
c  Several interfaces to the BMP He-H2 surface are provided.
c
c  Basic interface, as a function of interatomic distances:
c  ========================================================
c
c  call heh2bmp(r,Ra,Rb,id,V,dVdr,dVdRa,dVdRb,Vhh,dVhhdr,Vint,dVintdr)
c  -------------------------------------------------------------------
c
c  Input: r = H-H distance in bohrs
c         Ra, Rb = H-He distances in bohrs
c         id = (integer) derivative flag: id > 0 means "calculate the
c                analytic first derivatives", id = 0 means "do not do
c                so (return values of zero for all derivatives)"
c
c  Output: V = Vhh + Vint = total He-H2 surface energy in hartrees,
c                relative to separated H + H + He atoms
c          dVdr = derivative of V with respect to r at fixed Ra, Rb
c          dVdRa = derivative of V with respect to Ra at fixed r, Rb
c          dVdRb = derivative of V with respect to Rb at fixed r, Ra
c          Vhh = Schwenke H-H energy in hartrees as a function of r,
c                  relative to separated H + H atoms (note that this
c                  is modified at r < 0.1193892052031558 bohr, where
c                  Schwenke's formula fails, to yield 1/r at small r)
c          dVhhdr = derivative of Vhh with respect to r
c          Vint = remainder of He-H2 surface energy (if i_fix_close = 0
c                   then Vint = V_C + V_H + V_cb + V_d + V_d2; this may
c                   be modified at short distances if i_fix_close > 0:
c                   see description of heh2bmp_fix_close above)
c          dVintdr = derivative of Vint with respect to r at fixed Ra,
c                      Rb; note that dVintdRa and dVintdRb would be
c                      equal to dVdRa and dVdRb
c
c  Angular interface:
c  ==================
c
c  call heh2bmp_ang(r,RR,g,id,V,dr,da,db,Vh,dh,Vi,di,Ra,Rb,drgR,dRR,dg)
c  --------------------------------------------------------------------
c
c  Input: r = H-H distance
c         RR = H2 - He distance
c         g = H2 - He angle in radians
c         id = (integer) derivative flag
c
c  Output: V = Vh + Vi = total He-H2 surface energy
c          dr = derivative of V with respect to r at fixed Ra, Rb
c          da = derivative of V with respect to Ra at fixed r, Rb
c          db = derivative of V with respect to Rb at fixed r, Ra
c          Vh = Schwenke H-H energy as a function of r
c          dh = derivative of Vh with respect to r
c          Vi = remainder of He-H2 surface energy
c          di = derivative of Vi with respect to r at fixed Ra, Rb
c          Ra = sqrt( (RR*cos(g)-a2)**2 + (RR*sin(g))**2 )
c          Rb = sqrt( (RR*cos(g)+a2)**2 + (RR*sin(g))**2 )
c          drgR = derivative of V with respect to r at fixed RR and g
c          dRR = derivative of V with respect to RR at fixed r and g
c          dg = derivative of V with respect to g at fixed r and RR
c
c  General cartesian interface, including cartesian derivatives:
c  =============================================================
c
c  call heh2bmp_cc(c1,c2,ce,id,V,dr,da,db,Vh,dh,Vi,di,r,Ra,Rb,d1,d2,de)
c  --------------------------------------------------------------------
c
c  Input: c1(3) = cartesian coordinates (x1,y1,z1) of first H-atom
c         c2(3) = cartesian coordinates (x2,y2,z2) of second H-atom
c         ce(3) = cartesian coordinates (x3,y3,z3) of He-atom
c         id = (integer) derivative flag
c
c  Output: V = Vh + Vi = total He-H2 surface energy
c          dr = derivative of V with respect to r at fixed Ra, Rb
c          da = derivative of V with respect to Ra at fixed r, Rb
c          db = derivative of V with respect to Rb at fixed r, Ra
c          Vh = Schwenke H-H energy as a function of r
c          dh = derivative of Vh with respect to r
c          Vi = remainder of He-H2 surface energy
c          di = derivative of Vi with respect to r at fixed Ra, Rb
c          r = sqrt( (x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2 )
c          Ra = sqrt( (x3-x2)**2 + (y3-y2)**2 + (z3-z2)**2 )
c          Rb = sqrt( (x3-x1)**2 + (y3-y1)**2 + (z3-z1)**2 )
c          d1(3) = derivatives (dVdx1,dVdy1,dVdz1) of V with respect to
c                    cartesian coordinates of first H-atom
c          d2(3) = derivatives (dVdx2,dVdy2,dVdz2) of V with respect to
c                    cartesian coordinates of second H-atom
c          de(3) = derivatives (dVdx3,dVdy3,dVdz3) of V with respect to
c                    cartesian coordinates of He-atom
c
c  Note: dVdx1,...,dVdz3 are the derivatives with all the other eight
c        cartesian coordinates held fixed in each case.
c
c  Minimal cartesian interface (useful mainly for fitted geometries):
c  ==================================================================
c
c  call heh2bmp_cc3(a2,y,z,id,V,dr,da,db,Vh,dh,Vi,di,r,Ra,Rb,da2,dy,dz)
c  --------------------------------------------------------------------
c
c  Input: a2 = 0.5 * r = position on z-axis of the 2nd H-atom
c         y = y3 = distance of He-atom from z-axis
c         z = z3 = height of He-atom above x-y plane
c         id = (integer) derivative flag
c
c  Output: V = Vh + Vi = total He-H2 surface energy
c          dr = derivative of V with respect to r at fixed Ra, Rb
c          da = derivative of V with respect to Ra at fixed r, Rb
c          db = derivative of V with respect to Rb at fixed r, Ra
c          Vh = Schwenke H-H energy as a function of r
c          dh = derivative of Vh with respect to r
c          Vi = remainder of He-H2 surface energy
c          di = derivative of Vi with respect to r at fixed Ra, Rb
c          r = 2.0 * a2 = H-H distance
c          Ra = sqrt( (z3-a2)**2 + y3**2 )
c          Rb = sqrt( (z3+a2)**2 + y3**2 )
c          da2 = derivative of V with respect to a2 at fixed y and z
c          dy = derivative of V with respect to y at fixed a2 and z
c          dz = derivative of V with respect to z at fixed a2 and y
c
c   *** NOTE that da2, dy, and dz are not really meaningful.  Since
c          the positions of BOTH H-atoms depend on a2, these three
c          derivatives are NOT equal to the corresponding cartesian
c          coordinate derivatives in heh2bmp_cc above.
c
c######################################################################
c
c  Similar interfaces are supplied to the Muchnick and Russek surface:
c
c  call heh2MR(r,Ra,Rb,id,V,dVdr,dVdRa,dVdRb,Vhh,dVhhdr,Vint,dVintdr)
c  ------------------------------------------------------------------
c  call VintMR(r,Ra,Rb,id,Vint,dVintdr,dVdRa,dVdRb)
c  ------------------------------------------------
c
c  and to the "modMR" (modified Muchnick and Russek) He-H2 surface:
c
c  call heh2modMR(r,Ra,Rb,id,V,dVdr,dVdRa,dVdRb,Vhh,dVhhdr,Vint,dVintdr)
c  ---------------------------------------------------------------------
c  call VintmodMR(r,Ra,Rb,id,Vint,dVintdr,dVdRa,dVdRb)
c  ---------------------------------------------------
c
c  More general interfaces are supplied with an additional input "nD",
c  where nD = 0 yields the Muchnick and Russek surface and nD = 2 yields
c  the "modMR" surface (nD = 1 yields an intermediate surface, with V_d
c  as in the "modMR" surface but without the V_d2 = V_D term):
c
c  call heh2MR_nD(r,Ra,Rb,id,nD,V,dVdr,dVdRa,dVdRb,Vhh,dVhhdr,Vint,dVintdr)
c  ------------------------------------------------------------------------
c  call heh2MR_ang(r,RR,g,id,nD,V,dr,da,db,Vh,dh,Vi,di,Ra,Rb,drgR,dRR,dg)
c  ----------------------------------------------------------------------
c  call heh2MR_cc3(a2,y,z,id,nD,V,dr,da,db,Vh,dh,Vi,di,r,Ra,Rb,da2,dy,dz)
c  ----------------------------------------------------------------------
c  call heh2MR_cc(c1,c2,ce,id,nD,V,dr,da,db,Vh,dh,Vi,di,r,Ra,Rb,d1,d2,de)
c  ----------------------------------------------------------------------
c    with dimension c1(3), c2(3), ce(3), d1(3), d2(3), de(3)
c
c  and also interfaces to the Muchnick and Russek interaction portion:
c
c  call VintMR_nD(r,Ra,Rb,id,nD,Vint,dVintdr,dVdRa,dVdRb)
c  ------------------------------------------------------
c  call VintMR_ang(r,RR,g,id,nD,Vi,di,da,db,Ra,Rb,drgRi,dRR,dg)
c  ------------------------------------------------------------
c  call VintMR_cc3(a2,y,z,id,nD,Vi,di,da,db,r,Ra,Rb,da2i,dy,dz)
c  ------------------------------------------------------------
c  call VintMR_cc(c1,c2,ce,id,nD,Vi,di,da,db,r,Ra,Rb,d1i,d2i,de)
c  -------------------------------------------------------------
c    with dimension c1(3), c2(3), ce(3), d1i(3), d2i(3), de(3)
c
c######################################################################
c
c  An interface is supplied to the Schwenke H-H potential:
c
c  call Vhh_heh2bmp(r,id,Vhh,dVhhdr)
c  ---------------------------------
c
c  Input: r = H-H distance in bohrs
c         id = (integer) derivative flag: id > 0 means "calculate the
c              analytic first derivative", id = 0 means "do not do so"
c
c  Output: Vhh = Schwenke H-H potential energy in hartrees
c          dVhhdr = derivative of Vhh with respect to r (for id > 0)
c
c  Note that the actual Schwenke potential behaves poorly at small r,
c  so Vhh is switched over to 1/r, as r drops from 0.1193892052031558
c  to 0.0451564988896124 bohr; in consequence, the third derivative of
c  Vhh is discontinuous at r = 0.1193892052031558 bohr, and the second
c  derivative is discontinuous at r = 0.0451564988896124 bohr.
c
c
c**********************************************************************
c
c  Set the short-distance extrapolation flag:
c  ==========================================
c
      subroutine heh2bmp_fix_close( i_fix_close )
c     -------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
c__Note that a large input value of i_fix_close (9 or more) is reset to the
c  "best" switchover type (ifixclose = 2) for short-distance extrapolation.
c
      if ( i_fix_close .ge. 9 ) then
         ifixclose = 2
      else
         ifixclose = max( 0 , min( 3 , i_fix_close ) )
      endif
c
      return
      end
c
c**********************************************************************
c
c  Set the error level:
c  ====================
c
      subroutine heh2bmp_set_err( level )
c     -----------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
c__Note that large positive or negative inputs are reset to the (positive or
c  negative) default error level value.
c
      if ( level .ge. 9 ) then
         level_err = 1
      else if ( level .le. -9 ) then
         level_err = -1
      else
         level_err = max( -3 , min( 3 , level ) )
      endif
c
c__Set the minimum allowed distances according to the error level.
c
      if ( abs( level_err ) .le. 1 ) then
         Rabmin_err = 0.001d0
         rmin_err = 0.001d0
      else if ( abs( level_err ) .eq. 2 ) then
         Rabmin_err = 0.6d0
         rmin_err = 0.001d0
      else
         Rabmin_err = 0.6d0
         rmin_err = 0.4d0
      endif
c
      return
      end
c
c**********************************************************************
c
c  Return the current extrapolation and error settings:
c  ====================================================
c
      subroutine heh2bmp_ask_set(i_fix_close,level,Rab_min,r_min,V_err)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
      i_fix_close = ifixclose
      level = level_err
      Rab_min = Rabmin_err
      r_min = rmin_err
      V_err = Verror
c
      return
      end
c
c**********************************************************************
c
c  The BMP He-H2 surface as a function of the cartesian coordinates:
c  =================================================================
c
      subroutine heh2bmp_cc(ccH1,ccH2,ccHe,ideriv,V,dVdr,dVdRa,dVdRb,
     $         Vhh,dVhhdr,Vint,dVintdr,r,Ra,Rb,dVdccH1,dVdccH2,dVdccHe)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      dimension ccH1(3), ccH2(3), ccHe(3),
     $     dVdccH1(3), dVdccH2(3), dVdccHe(3)
c
c__Get the three distances as a function of the coordinates.
c
      r = sqrt( ( ccH1(1) - ccH2(1) )**2
     $     + ( ccH1(2) - ccH2(2) )**2 + ( ccH1(3) - ccH2(3) )**2 )
      Ra = sqrt( ( ccH2(1) - ccHe(1) )**2
     $     + ( ccH2(2) - ccHe(2) )**2 + ( ccH2(3) - ccHe(3) )**2 )
      Rb = sqrt( ( ccH1(1) - ccHe(1) )**2
     $     + ( ccH1(2) - ccHe(2) )**2 + ( ccH1(3) - ccHe(3) )**2 )
c
c__Get the BMP He-H2 surface as a function of these distances.
c
      call heh2bmp(r,Ra,Rb,ideriv,V,dVdr,dVdRa,dVdRb,
     $             Vhh,dVhhdr,Vint,dVintdr)
c
c__Get the derivatives as a function of the cartesian coordinates.
c
      if ( ideriv .gt. 0 ) then
c
         dVdccH1(1) = dVdr * ( ccH1(1) - ccH2(1) ) / r
     $        + dVdRb * ( ccH1(1) - ccHe(1) ) / Rb
         dVdccH1(2) = dVdr * ( ccH1(2) - ccH2(2) ) / r
     $        + dVdRb * ( ccH1(2) - ccHe(2) ) / Rb
         dVdccH1(3) = dVdr * ( ccH1(3) - ccH2(3) ) / r
     $        + dVdRb * ( ccH1(3) - ccHe(3) ) / Rb
         dVdccH2(1) = - dVdr * ( ccH1(1) - ccH2(1) ) / r
     $        + dVdRa * ( ccH2(1) - ccHe(1) ) / Ra
         dVdccH2(2) = - dVdr * ( ccH1(2) - ccH2(2) ) / r
     $        + dVdRa * ( ccH2(2) - ccHe(2) ) / Ra
         dVdccH2(3) = - dVdr * ( ccH1(3) - ccH2(3) ) / r
     $        + dVdRa * ( ccH2(3) - ccHe(3) ) / Ra
         dVdccHe(1) = - dVdRb * ( ccH1(1) - ccHe(1) ) / Rb
     $        - dVdRa * ( ccH2(1) - ccHe(1) ) / Ra
         dVdccHe(2) = - dVdRb * ( ccH1(2) - ccHe(2) ) / Rb
     $        - dVdRa * ( ccH2(2) - ccHe(2) ) / Ra
         dVdccHe(3) = - dVdRb * ( ccH1(3) - ccHe(3) ) / Rb
     $        - dVdRa * ( ccH2(3) - ccHe(3) ) / Ra
c
      else
c
         do i = 1, 3
            dVdccH1(i) = 0.d0
            dVdccH2(i) = 0.d0
            dVdccHe(i) = 0.d0
         enddo
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  The BMP He-H2 surface as a function of minimal cartesian coordinates:
c  =====================================================================
c
      subroutine heh2bmp_cc3(a2,y3,z3,ideriv,V,dVdr,dVdRa,dVdRb,Vhh,
     $                   dVhhdr,Vint,dVintdr,r,Ra,Rb,dVda2,dVdy3,dVdz3)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
c__Get the three distances.
c
      r = 2.d0 * a2
      Ra = sqrt( ( z3 - a2 )**2 + y3**2 )
      Rb = sqrt( ( z3 + a2 )**2 + y3**2 )
c
c__Get the BMP He-H2 surface as a function of these distances.
c
      call heh2bmp(r,Ra,Rb,ideriv,V,dVdr,dVdRa,dVdRb,
     $             Vhh,dVhhdr,Vint,dVintdr)
c
c__Get the derivatives as a function of the minimal cartesian coordinates
c  (these would not usually be useful).
c
      if ( ideriv .gt. 0 ) then
c
         dr_da2 = 2.d0
         dRa_da2 = - ( z3 - a2 ) / Ra
         dRb_da2 = ( z3 + a2 ) / Rb
         dRa_dy3 = y3 / Ra
         dRb_dy3 = y3 / Rb
         dRa_dz3 = ( z3 - a2 ) / Ra
         dRb_dz3 = ( z3 + a2 ) / Rb
c
         dVda2 = dVdr * dr_da2 + dVdRa * dRa_da2 + dVdRb * dRb_da2
         dVdy3 = dVdRa * dRa_dy3 + dVdRb * dRb_dy3
         dVdz3 = dVdRa * dRa_dz3 + dVdRb * dRb_dz3
c
      else
c
         dVda2 = 0.d0
         dVdy3 = 0.d0
         dVdz3 = 0.d0
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  The BMP He-H2 surface as a function of angular coordinates:
c  ===========================================================
c
      subroutine heh2bmp_ang(rhh,rh2he,gammarad,ideriv,V,dVdr,dVdRa,
     $     dVdRb,Vhh,dVhhdr,Vint,dVintdr,Ra,Rb,dVdrhh,dVdrh2he,dVdgam)
c     ----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_angderiv_BMP/ dRa_drhh, dRb_drhh,
     $     dRa_drh2he, dRb_drh2he, dRa_dgam, dRb_dgam
      save /vheh2fit_angderiv_BMP/
c
c__Get the distances Ra and Rb as a function of rhh, rh2he, and the
c  angle betwen them in radians gammarad.
c
      cgam = cos(gammarad)
      sgam = sin(gammarad)
      cgamR = cgam * rh2he
      sgamR = sgam * rh2he
      Ra = sqrt( ( cgamR - 0.5d0 * rhh )**2 + sgamR**2 )
      Rb = sqrt( ( cgamR + 0.5d0 * rhh )**2 + sgamR**2 )
c
c__Get the BMP He-H2 surface as a function of these distances.
c
      call heh2bmp(rhh,Ra,Rb,ideriv,V,dVdr,dVdRa,dVdRb,
     $             Vhh,dVhhdr,Vint,dVintdr)
c
c__Get the derivatives as a function of the angular coordinates.
c
      if ( ideriv .gt. 0 ) then
c
         dRa_drhh = -0.5d0 * ( cgamR - 0.5d0 * rhh ) / Ra
         dRb_drhh = 0.5d0 * ( cgamR + 0.5d0 * rhh ) / Rb
         dRa_drh2he = ( cgam * ( cgamR - 0.5d0 * rhh )
     $        + sgam * sgamR ) / Ra
         dRb_drh2he = ( cgam * ( cgamR + 0.5d0 * rhh )
     $        + sgam * sgamR ) / Rb
         dRa_dgam = ( cgamR * sgamR
     $        - sgamR * ( cgamR - 0.5d0 * rhh ) ) / Ra
         dRb_dgam = ( cgamR * sgamR
     $        - sgamR * ( cgamR + 0.5d0 * rhh ) ) / Rb
c
         dVdrhh = dVdr + dVdRa * dRa_drhh + dVdRb * dRb_drhh
         dVdrh2he = dVdRa * dRa_drh2he + dVdRb * dRb_drh2he
         dVdgam = dVdRa * dRa_dgam + dVdRb * dRb_dgam
c
      else
c
         dVdrhh = 0.d0
         dVdrh2he = 0.d0
         dVdgam = 0.d0
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  The BMP He-H2 surface as a function of the three distances:
c  ===========================================================
c
      subroutine heh2bmp( rhh,Ra,Rb, ideriv, V, dVdr,dVdRa,dVdRb,
     $                    Vhh, dVhhdr, Vint, dVintdr )
c     -----------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c					  ! upper limit for argument of exp() :
c-noneed;      parameter ( hi_exp=150.d0 )
c
c__Parameters used in long-range H-He interaction term V_d2:
c
      parameter ( harkcal=1.d0/627.5096d0, bohrangs=0.52917725d0 )
      parameter ( A4m_geng=-3.597516D-2, B4_geng=3.2D0 )
c
c__Largest possible index-limits in BMP-surface summations:
c
      parameter ( Lmax=6, Mmax=4, Imax=4, Imin=-Imax, Nmin=6, Nmax=12,
     $     Ncont=Nmin+2, Iminm1=Imin-1, Jmax=5,
     $     Jmaxp1=Jmax+1, Jmaxp2=Jmax+2 )
c
c__Actual index-limits in BMP-surface summations:
c
      common /vheh2fit_poly_BMP/ Mrac, Lac, Iac, Jac,
     $     MrH, LH, IH, JH, Mrcb, Lcb, Icb, Jcb, NRd, Ld, Mrd, NRheh,
     $     Mrmbe, Mhembe, MBEord
      save /vheh2fit_poly_BMP/
c
c__Fitted parameters of the BMP surface:
c
      parameter ( nparBMP = 500, ndelBMP = 22 )
c
      common /vheh2fit_par_BMP/ parBMP(nparBMP), NuseBMP
      save /vheh2fit_par_BMP/
      equivalence
     $     (parBMP(1),fac_A_d),    (parBMP(2),fac_alpha_d),
     $     (parBMP(3),e_ex),       (parBMP(4),d_ex),
     $     (parBMP(5),c_ex),       (parBMP(6),b_ex),
     $     (parBMP(7),c_heh),      (parBMP(8),C_1),
     $     (parBMP(9),C_2),        (parBMP(10),C_3),
     $     (parBMP(11),Acb_norm),  (parBMP(12),Acb_ex),
     $     (parBMP(13),dlcb_con),  (parBMP(14),dlcb_norm),
     $     (parBMP(15),dlcb_ex),   (parBMP(16),alpha_cb),
     $     (parBMP(17),A_H),       (parBMP(18),a),
     $     (parBMP(19),dlambda_p), (parBMP(20),AC_norm),
     $     (parBMP(21),AC_ex),     (parBMP(22),dlambda_c)
c
c__NOTE: fitted parameters parBMP(ipar) come in a number of independent sets:
c  set 0: basic BMP: ipar = 1, 22 :
c           fac_A_d, fac_alpha_d, e_ex, d_ex, c_ex, b_ex, c_heh, C_1, C_2, C_3,
c           Acb_norm, Acb_ex, dlcb_con, dlcb_norm, dlcb_ex, alpha_cb,
c           A_H, a, dlambda_p, AC_norm, AC_ex, dlambda_c
c  set 1: A_C poly: ( ( ( a_AC(i,L,m)=parBMP(ipar++), i = 0, Iac ),
c                           L = 0, Lac, 2 ), m = 1, Mrac )
c                   EXCEPT: skip case { i=0, L=0, m=1 }
c    ---for initial- or fit-parameter input purposes, set 1 is broken up into:
c            set 1: as above
c            set 11: beta1, pow1  for cutoff  exp( - beta1 * rho**pow1 )
c  set 2: V_C: polys of R_A, R_B: ( b_AC(j)=parBMP(ipar++), j = 1, Jac )
c  set 3: A_H poly: ( ( ( a_H(i,L,m)=parBMP(ipar++), i = 0, IH ),
c                           L = 0, LH, 2 ), m = 0, MrH )
c                   EXCEPT: skip case { i=0, L=0, m=0 }
c    ---for initial- or fit-parameter input purposes, set 3 is broken up into:
c            set 3: as above
c            set 12: beta3, pow3  for cutoff  exp( - beta3 * rho**pow3 )
c  set 4: V_H: polys of R_A, R_B: ( b_H(j)=parBMP(ipar++), j = 1, JH )
c  set 5: Acb poly: ( ( ( a_cb(i,L,m)=parBMP(ipar++), i = Icb, Jcb ),
c                           L = 0, Lcb, 2 ), m = 0, Mrcb )
c                   EXCEPT: skip cases { i=0, L=0, m=0 } and { i=0, L=2, m=0 }
c    ---for initial- or fit-parameter input purposes, set 5 is broken up into:
c            set 5: as above
c            set 13: beta5, pow5  for cutoff  exp( - beta5 * rho**pow5 )
c  set 6: V_d: ( ( a_d(m,L,6)=parBMP(ipar++), m = 1, Md(L,6) ),
c                     [ if ( Md(L,6) > 0 ) e_d(L,6)=parBMP(ipar++) ],
c                     L = 4, Ld, 2 ) ;
c              ( ( ( a_d(m,L,n)=parBMP(ipar++), m = min(L,1), Md(L,n) ),
c                     [ if ( Md(L,n) > min{L-1,0} ) e_d(L,n)=parBMP(ipar++) ],
c                     L = 0, Ld, 2 ),
c                     [ if ( max_over_L{ Md(L,n) - min(L-1,0) } > 0 )
c                          c_d(n)=parBMP(ipar++) ], n = 8, NRd, 2 )
c    ---for initial- or fit-parameter input purposes, set 6 is broken up into:
c            set 6: ( ( ( a_d(m,L,n), m = ... ), L = ... ), n = ... )
c            set 7: ( ( e_d(L,n), L = ... ), n = ... )
c            set 8: ( c_d(n), n = ... )
c  set 9: V_d(HeH): ( a_HeH(n)=parBMP(ipar++), c_HeH(n)=parBMP(ipar++),
c                       n = 8, NRheh, 2 )
c    ---for initial- or fit-parameter input purposes, set 9 is broken up into:
c            set 9: ( a_HeH(n), n = ... )
c            set 10: ( c_HeH(n), n = ... )
c
c__The terms in the BMP surface; note that, if ifixclose > 0, then at short
c  distances Vint_BMP may no longer be equal to V_C + V_H + V_cb + V_d + V_d2:
c
      common /vheh2fit_val_BMP/ Vhh_schwenke, Vint_BMP,
     $     V_C, V_P, V_d, V_d2, V_H, V_cb
      save /vheh2fit_val_BMP/
c
c__Derivatives of the terms in the BMP surface; note that dr_Vint0, dRa_Vint0,
c  and dRb_Vint0 are derivatives of Vint_0 = V_C + V_H + V_cb + V_d + V_d2,
c  not necessarily of Vint_BMP:
c
      common /vheh2fit_deriv_BMP/ d_Vhh,
     $     dr_Vint0, dRa_Vint0, dRb_Vint0, dr_V_C, dRa_V_C, dRb_V_C,
     $     dr_V_P, dRa_V_P, dRb_V_P, dr_V_d, dRa_V_d, dRb_V_d,
     $     dr_V_d2, dRa_V_d2, dRb_V_d2,
     $     dr_V_H, dRa_V_H, dRb_V_H, dr_V_cb, dRa_V_cb, dRb_V_cb
      save /vheh2fit_deriv_BMP/
c
c__Values and derivatives that may be used in the switchovers at short
c  distances for ifixclose > 0:
c
      common /vheh2fit_fV_fder_BMP/
     $     Vint_0, fnmr, dr_fnmr, omfnmr,
     $     r, R_A, R_B, fVint_Rab, ofVint_Rab, fcl2_Rab, Vcl, Vmr,
     $     dr_drhh, dr_drhhe1, dr_dRb,
     $     dR_A_drhh, dR_A_drhhe1, dR_A_dRb,
     $     dR_B_drhh, dR_B_dRa, dR_B_dRb,
     $     dfVint_drhh, dfVint_dRa, dfVint_dRb,
     $     dofVint_drhh, dofVint_dRa, dofVint_dRb,
     $     dfcl2_drhh, dfcl2_dRa, dfcl2_dRb,
     $     dr_Vcl, dR_A_Vcl, dR_B_Vcl,
     $     dr_Vmr, dRa_Vmr, dRb_Vmr
      save /vheh2fit_fV_fder_BMP/
c
c__The angular parameter  eta = ( Ra - Rb ) / r  and terms for the asymptotic
c  part of the long-range H2 - He interaction:
c
      common /vheh2fit_eta_BMP/ eta, alpha_d, alpha_d_eta
      save /vheh2fit_eta_BMP/
c
c__Legendre polynomials of eta, powers of the distances, and their derivatives:
c
      dimension PL_eta(0:Lmax), pow_Rbar(Iminm1:Nmax),
     $     pow_r(0:Jmaxp2), pow_Ra(-Jmaxp2:Nmax), pow_Rb(-Jmaxp2:Nmax),
     $     dr_PL_eta(0:Lmax), dRa_PL_eta(0:Lmax), dRb_PL_eta(0:Lmax),
     $     dRab_pow_Rbar(Imin:Nmax), d_pow_r(0:Jmaxp2),
     $     d_pow_Ra(-Jmaxp1:Nmax), d_pow_Rb(-Jmaxp1:Nmax)
      equivalence (PL_eta(2), P2_eta), (PL_eta(4), P4_eta),
     $     (PL_eta(6), P6_eta)
c
c__These would be used to store quantities for the (unused) MBE term:
c
c-noneed;      dimension pow_pr(0:Lmax), pow_pa(0:Lmax), pow_pb(0:Lmax)
c
c__Constants needed for computation of the Legendre polynomials:
c
      parameter ( P0_0 = 1.d0, P2_2 = 1.5d0, P2_0 = -0.5d0,
     $     P4_4 = 4.375d0, P4_2 = -3.75d0, P4_0 = 0.375d0,
     $     P6_6 = 14.4375d0, P6_4 = -19.6875d0,
     $     P6_2 = 6.5625d0, P6_0 = -0.3125d0 )
c
c__Constants needed for computation of long-range factors A_d and alpha_d:
c
      parameter ( f2_alpha_d = 0.238248d0, f1_alpha_d = -0.126241d0,
     $     fex_alpha_d = -0.8075d0, f2_A_d = 4.790521d0,
     $     f1_A_d = -2.494444d0, f0_A_d = 3.008850d0,
     $     fex_A_d = -0.5891d0 )
c
c__Values used to control error-handling and short-distance extrapolation:
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
c__For R_HeH = 0.5, VHeH_abinitio = 1.282424 = 2 * exp(heh_lam*0.5) / 0.5
c  (this is not used for the BMP surface):
c
c-noneed;      parameter ( heh_lam = -2.27508465d0 )
c
c__For R_HeH = 0.6, VHeH_abinitio = 0.874802 = 2 * exp(heh_cl2*0.6) / 0.6 ;
c  used for ifixclose = 2 or 3 in the factor used to estimate the He-H energy
c  at short distances Ra or Rb in terms of the energy at (or near) 0.6 bohr:
c
      parameter ( heh_cl2 = -2.22955085d0 )
c
c__For R_HeH = 0.7, VHeH_abinitio = 0.628941 = 2 * exp(heh_cl*0.7) / 0.7 ,
c  used for ifixclose = 1 to estimate the He-H energy at short distances Ra or
c  Rb when the value of r is relatively large (instead of the "modMR" energy);
c  the other  ..._cl  values determine the switchover point where Ra,Rb are
c  considered "short", and the  ..._mr  values determine the switchover from
c  the "modMR" energy as a function of r: 
c
      parameter ( heh_cl = -2.16219993d0,
     $     rlo_cl = 0.6d0, rhi_cl = 0.7d0, rmid_cl = 0.65d0,
     $     dinv_cl = 20.d0, oinv_cl = 10.d0,
     $     rlo_mr = 0.d0, rhi_mr = 8.d0, rmid_mr = 4.d0,
     $     dinv_mr = 0.25d0, oinv_mr = 0.125d0 )
c
c__For ifixclose = 2 or 3, these are the factors determining the switchover
c  to the alternate extrapolation at short Ra or Rb:
c
      parameter ( rmid_cl2 = 0.6d0, d_cl2 = 0.05d0,
     $     rlo_cl2 = rmid_cl2 - d_cl2, rhi_cl2 = rmid_cl2 + d_cl2,
     $     denom_cl2 = 4.d0 * d_cl2 * rlo_cl2 * rlo_cl2,
     $     b_cl2 = -1.d0 / denom_cl2,
     $     a_cl2 = ( rmid_cl2 - 3.d0 * d_cl2 ) / denom_cl2 )
c
c__For ifixclose = 3, these are the factors controlling the switchover to the
c  "modMR" surface at small r:
c
      parameter ( rlo_clh = 0.3d0, rhi_clh = 0.5d0, rmid_clh = 0.4d0,
     $     dinv_clh = 10.d0, oinv_clh = 5.d0 )
c
c-debug[
c-debug;      common /vheh2fit_debug/ idebug, ndebug, modebug
c-debug;      save /vheh2fit_debug/
c-debug]
c
c__Initialize all the derivatives to zero:
c
      dVhhdr = 0.d0
      dVintdr = 0.d0
      dVdr = 0.d0
      dVdRa = 0.d0
      dVdRb = 0.d0
c
c__Get local variables holding the three interatomic distance values (note that
c  these may be modified later, if ifixclose = 2 or 3):
c
      r = rhh
      R_A = Ra
      R_B = Rb
c
c__Check for validity of the geometry:
c
      call geomchk_heh2bmp( r, R_A, R_B, V, ierr )
c
c__If it is invalid, set all energies to the "error-value" 1.0E+36 and return:
c
      if ( ierr .ne. 0 ) then
         Vhh = V
         Vint = V
         return
      endif
c
c__Get the Schwenke H-H interaction energy:
c
      call Vhh_heh2bmp( r, ideriv, Vhh_schwenke, d_Vhh )
c
      Vhh = Vhh_schwenke
c
c__Initialize switchover values (may be modified later if ifixclose > 0):
c
      Vcl = 0.d0
      dr_Vcl = 0.d0
      dR_A_Vcl = 0.d0
      dR_B_Vcl = 0.d0
c
      Vmr = 0.d0
      dr_Vmr = 0.d0
      dRa_Vmr = 0.d0
      dRb_Vmr = 0.d0
c
      fnmr = 1.d0
      omfnmr = 0.d0
      dr_fnmr = 0.d0
c
      fVint_Rab = 1.d0
      dfVint_drhh = 0.d0
      dfVint_dRa = 0.d0
      dfVint_dRb = 0.d0
c
      ofVint_Rab = 0.d0
      dofVint_drhh = 0.d0
      dofVint_dRa = 0.d0
      dofVint_dRb = 0.d0
c
      dr_drhh = 1.d0
      dr_dRa = 0.d0
      dr_dRb = 0.d0
      dR_A_drhh = 0.d0
      dR_A_dRa = 1.d0
      dR_A_dRb = 0.d0
      dR_B_drhh = 0.d0
      dR_B_dRa = 0.d0
      dR_B_dRb = 1.d0
c
      r_lo = min( R_A , R_B )
c
c__If ifixclose = 1 and the shortest H-He distance is small enough that the
c  short-distance extrapolation should be used, compute relevant quantities:
c
      if ( ifixclose .eq. 1 .and. r_lo .lt. rhi_cl ) then
c							  ! if need only modMR:
         if ( r .le. rlo_mr ) then
c					! set switch factor and get Vint_modMR
            omfnmr = 1.d0
c-noneed;            dr_fnmr = 0.d0
c
            call VintmodMR(r,R_A,R_B,ideriv,Vmr,dr_Vmr,dRa_Vmr,dRb_Vmr)
c									! else:
         else
c		! computed expontially damped 2/R factors for V_cl
c
            Vcl_A = 2.d0 * exp( heh_cl * R_A ) / R_A
            Vcl_B = 2.d0 * exp( heh_cl * R_B ) / R_B
            Vcl = Vcl_A + Vcl_B
c					! and their derivatives
            if ( ideriv .gt. 0 ) then
               dR_A_Vcl = ( heh_cl - 1.d0 / R_A ) * Vcl_A
               dR_B_Vcl = ( heh_cl - 1.d0 / R_B ) * Vcl_B
            endif
c					! if modMR is needed too:
            if ( r .lt. rhi_mr ) then
c						! get Vint_modMR
               call VintmodMR(r,R_A,R_B,ideriv,
     $              Vmr,dr_Vmr,dRa_Vmr,dRb_Vmr)
c						! get modMR-switch factor
               if ( r .le. rmid_mr ) then
                  o_mxd = 0.5d0 - ( r - rmid_mr ) * oinv_mr
                  o_pxd = 1.d0 + ( r - rmid_mr ) * dinv_mr
                  fnmr = o_pxd**3 * o_mxd
                  omfnmr = 1.d0 - fnmr
                  if ( ideriv .gt. 0 ) dr_fnmr = oinv_mr
     $                 * o_pxd**2 * ( 6.d0 * o_mxd - o_pxd )
               else
                  o_mxd = 1.d0 - ( r - rmid_mr ) * dinv_mr
                  o_pxd = 0.5d0 + ( r - rmid_mr ) * oinv_mr
                  omfnmr = o_mxd**3 * o_pxd
                  fnmr = 1.d0 - omfnmr
                  if ( ideriv .gt. 0 ) dr_fnmr = oinv_mr
     $                 * o_mxd**2 * ( 6.d0 * o_pxd - o_mxd )
               endif
c
            endif
c
         endif
c					! get short-distance switch factor:
         if ( r_lo .le. rlo_cl ) then
c					! fully switched over
            fVint_Rab = 0.d0
            ofVint_Rab = 1.d0
c					! or only partially:
         else
c						! factor for Ra:
            if ( R_A .le. rmid_cl ) then
               o_mxd = 0.5d0 - ( R_A - rmid_cl ) * oinv_cl
               o_pxd = 1.d0 + ( R_A - rmid_cl ) * dinv_cl
               fcl_A = o_pxd**3 * o_mxd
               if ( ideriv .gt. 0 ) d_fcl_A = oinv_cl
     $              * o_pxd**2 * ( 6.d0 * o_mxd - o_pxd )
            else if ( R_A .lt. rhi_cl ) then
               o_mxd = 1.d0 - ( R_A - rmid_cl ) * dinv_cl
               o_pxd = 0.5d0 + ( R_A - rmid_cl ) * oinv_cl
               fcl_A = 1.d0 - o_mxd**3 * o_pxd
               if ( ideriv .gt. 0 ) d_fcl_A = oinv_cl
     $              * o_mxd**2 * ( 6.d0 * o_pxd - o_mxd )
            else
               fcl_A = 1.d0
               d_fcl_A = 0.d0
            endif
c						! factor for Rb:
            if ( R_B .le. rmid_cl ) then
               o_mxd = 0.5d0 - ( R_B - rmid_cl ) * oinv_cl
               o_pxd = 1.d0 + ( R_B - rmid_cl ) * dinv_cl
               fcl_B = o_pxd**3 * o_mxd
               if ( ideriv .gt. 0 ) d_fcl_B = oinv_cl
     $              * o_pxd**2 * ( 6.d0 * o_mxd - o_pxd )
            else if ( R_B .lt. rhi_cl ) then
               o_mxd = 1.d0 - ( R_B - rmid_cl ) * dinv_cl
               o_pxd = 0.5d0 + ( R_B - rmid_cl ) * oinv_cl
               fcl_B = 1.d0 - o_mxd**3 * o_pxd
               if ( ideriv .gt. 0 ) d_fcl_B = oinv_cl
     $              * o_mxd**2 * ( 6.d0 * o_pxd - o_mxd )
            else
               fcl_B = 1.d0
               d_fcl_B = 0.d0
            endif
c						! combine the two:
            fVint_Rab = fcl_A * fcl_B
            ofVint_Rab = 1.d0 - fVint_Rab
c
            if ( ideriv .gt. 0 ) then
c-noneed;               dfVint_drhh = 0.d0
               dfVint_dRa = d_fcl_A * fcl_B
               dfVint_dRb = fcl_A * d_fcl_B
c-noneed;               dofVint_drhh = 0.d0
               dofVint_dRa = -dfVint_dRa
               dofVint_dRb = -dfVint_dRb
            endif
c
         endif
c
c__ELSE: if ifixclose = 2 or 3 and the shortest H-He distance is small enough
c  that short-distance extrapolation is needed, compute relevant quantities:
c
      else if ( ifixclose .ge. 2 .and. r_lo .lt. rhi_cl2 ) then
c
c					! short-distance switch factor for Ra:
         if ( R_A .le. rlo_cl2 ) then
            fcl_A = rmid_cl2 / R_A - 1.d0
            d_fcl_A = - rmid_cl2 / R_A**2
         else if ( R_A .lt. rhi_cl2 ) then
            fcl_A = ( ( R_A - rhi_cl2 ) * b_cl2 + a_cl2 )
     $           * ( R_A - rhi_cl2 )**2
            d_fcl_A = ( 3.d0 * ( R_A - rhi_cl2 ) * b_cl2
     $           + 2.d0 * a_cl2 ) * ( R_A - rhi_cl2 )
         else
            fcl_A = 0.d0
            d_fcl_A = 0.d0
         endif
c					! short-distance switch factor for Rb:
         if ( R_B .le. rlo_cl2 ) then
            fcl_B = rmid_cl2 / R_B - 1.d0
            if ( ideriv .gt. 0 ) d_fcl_B = - rmid_cl2 / R_B**2
         else if ( R_B .lt. rhi_cl2 ) then
            fcl_B = ( ( R_B - rhi_cl2 ) * b_cl2 + a_cl2 )
     $           * ( R_B - rhi_cl2 )**2
            if ( ideriv .gt. 0 ) d_fcl_B =
     $           ( 3.d0 * ( R_B - rhi_cl2 ) * b_cl2
     $           + 2.d0 * a_cl2 ) * ( R_B - rhi_cl2 )
         else
            fcl_B = 0.d0
            d_fcl_B = 0.d0
         endif
c					! combine the two:
         if ( fcl_B .eq. 0.d0 ) then
            fcl2m1 = fcl_A
            dfcl2_dRa = d_fcl_A
            dfcl2_dRb = 0.d0
         else if ( fcl_A .eq. 0.d0 ) then
            fcl2m1 = fcl_B
            dfcl2_dRa = 0.d0
            dfcl2_dRb = d_fcl_B
         else
            fcl2m1 = ( fcl_A**10 + fcl_B**10 )**0.1d0
            if ( ideriv .gt. 0 ) then
               dfcl2_dRa = ( fcl_A / fcl2m1 )**9 * d_fcl_A
               dfcl2_dRb = ( fcl_B / fcl2m1 )**9 * d_fcl_B
            endif
         endif
c					! multiplicative factor for distances:
         fcl2_Rab = 1.d0 + fcl2m1
c					! get expanded distances: min{R_A,R_B}
c					! will not exceed 0.6 bohr now:
         r = rhh * fcl2_Rab
         R_A = Ra * fcl2_Rab
         R_B = Rb * fcl2_Rab
c					    ! multiplicative factor for Vint:
         den_A = exp( heh_cl2 * R_A ) / R_A
         den_B = exp( heh_cl2 * R_B ) / R_B 
         denom = den_A + den_B
         Vnum_1 = exp( heh_cl2 * Ra ) / Ra
         Vnum_2 = exp( heh_cl2 * Rb ) / Rb
         fVint_Rab = ( Vnum_1 + Vnum_2 ) / denom
c
c-noneed;         ofVint_Rab = 0.d0
c					! compute the necessary derivatives:
         if ( ideriv .gt. 0 ) then
c
c-noneed;            dfcl2_drhh = 0.d0
c					! of new distances relative to the old:
            dr_drhh = fcl2_Rab
            dr_dRa = rhh * dfcl2_dRa
            dr_dRb = rhh * dfcl2_dRb
c-noneed;            dR_A_drhh = 0.d0
            dR_A_dRa = fcl2_Rab + Ra * dfcl2_dRa
            dR_A_dRb = Ra * dfcl2_dRb
c-noneed;            dR_B_drhh = 0.d0
            dR_B_dRa = Rb * dfcl2_dRa
            dR_B_dRb = fcl2_Rab + Rb * dfcl2_dRb
c						  ! and of factor for Vint:
c-noneed;            dfVint_drhh = 0.d0
            dfVint_dRa =
     $           ( ( 1.d0 / R_A - heh_cl2 ) * den_A * dR_A_dRa
     $           + ( 1.d0 / R_B - heh_cl2 ) * den_B * dR_B_dRa )
     $           * fVint_Rab / denom
     $           + ( heh_cl2 - 1.d0 / Ra ) * Vnum_1 / denom
            dfVint_dRb =
     $           ( ( 1.d0 / R_A - heh_cl2 ) * den_A * dR_A_dRb
     $           + ( 1.d0 / R_B - heh_cl2 ) * den_B * dR_B_dRb )
     $           * fVint_Rab / denom
     $           + ( heh_cl2 - 1.d0 / Rb ) * Vnum_2 / denom
c
c-noneed;            dofVint_drhh = 0.d0
c-noneed;            dofVint_dRa = 0.d0
c-noneed;            dofVint_dRb = 0.d0
c
         endif
c
      endif
c
c__If ifixclose = 3 and the H-H distance r is small enough that the "modMR"
c  switchover should be used, compute relevant quantities:
c
      if ( ifixclose .ge. 3 .and. rhh .lt. rhi_clh ) then
c							   ! swtich on modMR:
         fnmr = 0.d0
         omfnmr = 1.d0
c					 ! get Vint_modMR of original distances
         call VintmodMR(rhh,Ra,Rb,ideriv,
     $        Vmr,dr_Vmr,dRa_Vmr,dRb_Vmr)
c					  ! switch factor (vs. H-H distance):
         if ( rhh .le. rlo_clh ) then
            f_mul = 0.d0
            ofVint_Rab = 1.d0
c-noneed;            dofVint_drhh = 0.d0
         else if ( rhh .le. rmid_clh ) then
            o_mxd = 0.5d0 - ( rhh - rmid_clh ) * oinv_clh
            o_pxd = 1.d0 + ( rhh - rmid_clh ) * dinv_clh
            f_mul = o_pxd**3 * o_mxd
            ofVint_Rab = 1.d0 - f_mul
            if ( ideriv .gt. 0 ) dofVint_drhh = oinv_clh
     $           * o_pxd**2 * ( o_pxd - 6.d0 * o_mxd )
         else
            o_mxd = 1.d0 - ( rhh - rmid_clh ) * dinv_clh
            o_pxd = 0.5d0 + ( rhh - rmid_clh ) * oinv_clh
            ofVint_Rab = o_mxd**3 * o_pxd
            f_mul = 1.d0 - ofVint_Rab
            if ( ideriv .gt. 0 ) dofVint_drhh = oinv_clh
     $           * o_mxd**2 * ( o_mxd - 6.d0 * o_pxd )
         endif
c					! apply switch factor to that for Vint:
         if ( ideriv .gt. 0 ) then
            dfVint_drhh = dfVint_drhh * f_mul
     $           - fVint_Rab * dofVint_drhh
            dfVint_dRa = dfVint_dRa * f_mul
            dfVint_dRb = dfVint_dRb * f_mul
         endif
c
         fVint_Rab = fVint_Rab * f_mul
c
      endif
c
c__Get the spheroidal coordinates eta, Rbar used in the surface equations:
c
      eta = ( R_A - R_B ) / r
      Rbar = ( R_A + R_B ) * 0.5d0
c						! (exponential cutoff unused):
c-noneed;      rhosq = r**2 + R_A**2 + R_B**2
c
c-debug[
c-debug;      if ( idebug .gt. 0 .and. mod(idebug,modebug) .eq. 0 ) then
c-debug;         write(6,
c-debug;     $        '(" "/i5,": r,R_A,B",3f15.9,"  Rbar,eta",2f15.9,":")')
c-debug;     $        idebug,r, R_A, R_B, Rbar, eta
c-debug;      endif
c-debug]
c
c__Calculate Legendre functions of eta:
c
      etasq = eta**2
      PL_eta(0) = P0_0
      P2_eta = P2_2 * etasq + P2_0
      P4_eta = ( P4_4 * etasq + P4_2 ) * etasq + P4_0
      P6_eta = ( ( P6_6 * etasq + P6_4 ) * etasq
     $     + P6_2 ) * etasq + P6_0
c
c__Calculate the required powers of the distances r, R_A, R_B, Rbar:
c
      pow_r(0) = 1.d0
      pow_Ra(0) = 1.d0
      pow_Rb(0) = 1.d0
      pow_Rbar(0) = 1.d0
c
      do i = 1, Imax
         pow_Rbar(i) = pow_Rbar(i-1) * Rbar
         pow_Rbar(-i) = pow_Rbar(1-i) / Rbar
         pow_Ra(i) = pow_Ra(i-1) * R_A
         pow_Rb(i) = pow_Rb(i-1) * R_B
      enddo
c
c__IF derivatives are not to be calculated, fewer powers are needed:
c
      if ( ideriv .le. 0 ) then
c
         do i = Imax + 2, Nmax, 2
            pow_Rbar(i) = pow_Rbar(i-2) * pow_Rbar(2)
            pow_Ra(i) = pow_Ra(i-2) * pow_Ra(2)
            pow_Rb(i) = pow_Rb(i-2) * pow_Rb(2)
         enddo
c				! note Ra,b^-m are defined down to -6
         do m = 1, Jmaxp1
            pow_r(m) = pow_r(m-1) * r
            pow_Ra(-m) = pow_Ra(1-m) / R_A
            pow_Rb(-m) = pow_Rb(1-m) / R_B
         enddo
c
c__ELSE: derivatives are to be calculated: more powers are needed, and also
c  derivatives of the Legendre functions and of the powers of distances:
c
      else
c
         dr_etasq = -2.d0 * etasq / r
         dRa_etasq = 2.d0 * eta / r
         dRb_etasq = - dRa_etasq
c				 ! d PL / d r ,  d PL / d Ra ,  d PL / d Rb
         dr_PL_eta(0) = 0.d0
         dRa_PL_eta(0) = 0.d0
         dRb_PL_eta(0) = 0.d0
         dr_PL_eta(2) = P2_2 * dr_etasq
         dRa_PL_eta(2) = P2_2 * dRa_etasq
         dRb_PL_eta(2) = P2_2 * dRb_etasq
         dPL = 2.d0 * P4_4 * etasq + P4_2
         dr_PL_eta(4) = dPL * dr_etasq
         dRa_PL_eta(4) = dPL * dRa_etasq
         dRb_PL_eta(4) = dPL * dRb_etasq
         dPL = ( 3.d0 * P6_6 * etasq + 2.d0 * P6_4 ) * etasq + P6_2
         dr_PL_eta(6) = dPL * dr_etasq
         dRa_PL_eta(6) = dPL * dRa_etasq
         dRb_PL_eta(6) = dPL * dRb_etasq
c
         do i = Imax + 1, Nmax
            pow_Rbar(i) = pow_Rbar(i-1) * Rbar
            pow_Ra(i) = pow_Ra(i-1) * R_A
            pow_Rb(i) = pow_Rb(i-1) * R_B
         enddo
c
         pow_Rbar(Iminm1) = pow_Rbar(Imin) / Rbar
c							! d Rbar^i / d Ra,b
         do i = Imin, 0
            dRab_pow_Rbar(i) = 0.5d0 * i * pow_Rbar(i-1)
         enddo
         do i = 1, Nmax
            dRab_pow_Rbar(i) = 0.5d0 * i * pow_Rbar(i-1)
            d_pow_Ra(i) = i * pow_Ra(i-1)
            d_pow_Rb(i) = i * pow_Rb(i-1)
         enddo
c
         d_pow_r(0) = 0.d0
c				! note d Ra,b^-m down to -6 (but Ra,b^-m to -7)
         do m = 1, Jmaxp2
            pow_r(m) = pow_r(m-1) * r
            d_pow_r(m) = m * pow_r(m-1)
            pow_Ra(-m) = pow_Ra(1-m) / R_A
            pow_Rb(-m) = pow_Rb(1-m) / R_B
            d_pow_Ra(1-m) = ( 1 - m ) * pow_Ra(-m)
            d_pow_Rb(1-m) = ( 1 - m ) * pow_Rb(-m)
         enddo
c
      endif
c
c__Start the fitted-surface-parameter index  ipar  at the end of the basic
c  parameters, in preparation for calculating the V_C term:
c
      ipar = ndelBMP
c
c__Get the summation required for A_C (used in the V_C term):
c
      A_C_poly = 0.d0
      jpar = ipar
      ilo = -1
c
      if ( ideriv .le. 0 ) then
c
         do m = 1, Mrac
            vaddm = 0.d0
            do L = 0, Lac, 2
               vadd = 0.d0
               do i = ilo, -Iac, -1
                  ipar = ipar + 1
                  vadd = vadd + parBMP(ipar) * pow_Rbar(i)
               enddo
               vaddm = vaddm + PL_eta(L) * vadd
               ilo = 0
            enddo
            A_C_poly = A_C_poly + pow_r(m) * vaddm
         enddo
c
      else
c
         dr_A_C_poly = 0.d0
         dRa_A_C_poly = 0.d0
         dRb_A_C_poly = 0.d0
c
         do m = 1, Mrac
            vaddm = 0.d0
            dr_vaddm = 0.d0
            dRa_vaddm = 0.d0
            dRb_vaddm = 0.d0
            do L = 0, Lac, 2
               vadd = 0.d0
               d_vadd = 0.d0
               do i = ilo, -Iac, -1
                  ipar = ipar + 1
                  vadd = vadd + parBMP(ipar) * pow_Rbar(i)
                  d_vadd = d_vadd + parBMP(ipar) * dRab_pow_Rbar(i)
               enddo
               vaddm = vaddm + PL_eta(L) * vadd
               dr_vaddm = dr_vaddm + dr_PL_eta(L) * vadd
               dRa_vaddm = dRa_vaddm + dRa_PL_eta(L) * vadd
     $              + PL_eta(L) * d_vadd
               dRb_vaddm = dRb_vaddm + dRb_PL_eta(L) * vadd
     $              + PL_eta(L) * d_vadd
               ilo = 0
            enddo
            A_C_poly = A_C_poly + pow_r(m) * vaddm
            dr_A_C_poly = dr_A_C_poly + d_pow_r(m) * vaddm
     $           + pow_r(m) * dr_vaddm
            dRa_A_C_poly = dRa_A_C_poly + pow_r(m) * dRa_vaddm
            dRb_A_C_poly = dRb_A_C_poly + pow_r(m) * dRb_vaddm
         enddo
c
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'aAC: ipar=',ipar
c-debug]
c
c__Skip over the parameters controlling the unused exponential cutoff:
c
      if ( ipar .gt. jpar ) then
         ipar = ipar + 1
         beta = parBMP(ipar)
         ipar = ipar + 1
         if ( beta .ne. 0.d0 ) then
            stop ' STOP -- Error: exp(beta*rho) term not implemented. '
c-noneed;            pow = parBMP(ipar)
c-noneed;            if ( pow .eq. 0.d0 ) then
c-noneed;               rhop = 1.d0
c-noneed;            else if ( pow .eq. 2.d0 ) then
c-noneed;               rhop = rhosq
c-noneed;            else
c-noneed;               rhop = rhosq**( 0.5d0 * pow )
c-noneed;            endif
c-noneed;            A_C_poly = A_C_poly * exp( min( - beta * rhop , hi_exp ) )
         endif
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'BAC: ipar=',ipar
c-debug]
c
c__Compute the final value of A_C (for the V_C term):
c
c-noneed;      A_C = 2.d0 - ( AC_norm * r + A_C_poly )
c-noneed;     $     * exp( min( AC_ex * r , hi_exp ) )
c
      A_C_exfac = exp( AC_ex * r )
      A_C = 2.d0 - ( AC_norm * r + A_C_poly ) * A_C_exfac
c
      if ( ideriv .gt. 0 ) then
c
         dr_A_C = - AC_ex * ( AC_norm * r + A_C_poly ) * A_C_exfac
     $        - ( AC_norm + dr_A_C_poly ) * A_C_exfac
         dRa_A_C = - dRa_A_C_poly * A_C_exfac
         dRb_A_C = - dRb_A_C_poly * A_C_exfac
c
      endif
c
c__Compute the F_C summations that depend on R_A or R_B alone (for V_C term):
c
      R_A_poly = pow_Ra(-1)
      R_B_poly = pow_Rb(-1)
c
      if ( ideriv .le. 0 ) then
c
         do j = -2, -1 - Jac, -1
            ipar = ipar + 1
            R_A_poly = R_A_poly + parBMP(ipar) * pow_Ra(j)
            R_B_poly = R_B_poly + parBMP(ipar) * pow_Rb(j)
         enddo
c
      else
c
         d_R_A_poly = d_pow_Ra(-1)
         d_R_B_poly = d_pow_Rb(-1)
c
         do j = -2, -1 - Jac, -1
            ipar = ipar + 1
            R_A_poly = R_A_poly + parBMP(ipar) * pow_Ra(j)
            R_B_poly = R_B_poly + parBMP(ipar) * pow_Rb(j)
            d_R_A_poly = d_R_A_poly + parBMP(ipar) * d_pow_Ra(j)
            d_R_B_poly = d_R_B_poly + parBMP(ipar) * d_pow_Rb(j)
         enddo
c
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'bAC: ipar=',ipar
c-debug]
c
c__Compute the V_C term and its derivatives, using A_C and the F_C summations:
c
c-noneed;      V_C = A_C * ( 1.d0 + ( 2.d0 / A_C - 1.d0 ) * P2_eta )
c-noneed;     $     * ( exp( min( - dlambda_c * R_A , hi_exp ) ) * R_A_poly
c-noneed;     $     + exp( min( - dlambda_c * R_B , hi_exp ) ) * R_B_poly )
c
      f_C_Ra = exp( - dlambda_c * R_A ) * R_A_poly
      f_C_Rb = exp( - dlambda_c * R_B ) * R_B_poly
c
      A_C_ang = 1.d0 + ( 2.d0 / A_C - 1.d0 ) * P2_eta
      V_C = A_C * A_C_ang * ( f_C_Ra + f_C_Rb )
c
      if ( ideriv .gt. 0 ) then
c
         dr_V_C = ( dr_A_C * A_C_ang
     $        + A_C * ( ( 2.d0 / A_C - 1.d0 ) * dr_PL_eta(2)
     $        - 2.d0 * dr_A_C / A_C**2 * P2_eta ) )
     $        * ( f_C_Ra + f_C_Rb )
         dRa_V_C = ( dRa_A_C * A_C_ang
     $        + A_C * ( ( 2.d0 / A_C - 1.d0 ) * dRa_PL_eta(2)
     $        - 2.d0 * dRa_A_C / A_C**2 * P2_eta ) )
     $        * ( f_C_Ra + f_C_Rb )
     $        + A_C * A_C_ang * ( exp( - dlambda_c * R_A ) * d_R_A_poly
     $        - dlambda_c * f_C_Ra )
         dRb_V_C = ( dRb_A_C * A_C_ang
     $        + A_C * ( ( 2.d0 / A_C - 1.d0 ) * dRb_PL_eta(2)
     $        - 2.d0 * dRb_A_C / A_C**2 * P2_eta ) )
     $        * ( f_C_Ra + f_C_Rb )
     $        + A_C * A_C_ang * ( exp( - dlambda_c * R_B ) * d_R_B_poly
     $        - dlambda_c * f_C_Rb )
c
      endif
c
c-debug[
c-debug;      if ( idebug .gt. 0 .and. mod(idebug,modebug) .eq. 0 ) then
c-debug;         write(6,'(" ipar",i4,"  A_C",f15.9,"  V_C",f15.9)')
c-debug;     $        ipar, A_C, V_C
c-debug;      endif
c-debug]
c
c__Compute the summation for the V_H term:
c
      a_H_poly = 0.d0
c
      jpar = ipar
      ilo = -1
c
      if ( ideriv .le. 0 ) then
c
         do m = 0, MrH
            vaddm = 0.d0
            do L = 0, LH, 2
               vadd = 0.d0
               do i = ilo, -IH, -1
                  ipar = ipar + 1
                  vadd = vadd + parBMP(ipar) * pow_Rbar(i)
               enddo
               vaddm = vaddm + PL_eta(L) * vadd
               ilo = 0
            enddo
            a_H_poly = a_H_poly + pow_r(m) * vaddm
         enddo
c
      else
c
         dr_a_H_poly = 0.d0
         dRa_a_H_poly = 0.d0
         dRb_a_H_poly = 0.d0
c
         do m = 0, MrH
            vaddm = 0.d0
            dr_vaddm = 0.d0
            dRa_vaddm = 0.d0
            dRb_vaddm = 0.d0
            do L = 0, LH, 2
               vadd = 0.d0
               d_vadd = 0.d0
               do i = ilo, -IH, -1
                  ipar = ipar + 1
                  vadd = vadd + parBMP(ipar) * pow_Rbar(i)
                  d_vadd = d_vadd + parBMP(ipar) * dRab_pow_Rbar(i)
               enddo
               vaddm = vaddm + PL_eta(L) * vadd
               dr_vaddm = dr_vaddm + dr_PL_eta(L) * vadd
               dRa_vaddm = dRa_vaddm + dRa_PL_eta(L) * vadd
     $              + PL_eta(L) * d_vadd
               dRb_vaddm = dRb_vaddm + dRb_PL_eta(L) * vadd
     $              + PL_eta(L) * d_vadd
               ilo = 0
            enddo
            a_H_poly = a_H_poly + pow_r(m) * vaddm
            dr_a_H_poly = dr_a_H_poly + pow_r(m) * dr_vaddm
     $           + d_pow_r(m) * vaddm
            dRa_a_H_poly = dRa_a_H_poly + pow_r(m) * dRa_vaddm
            dRb_a_H_poly = dRb_a_H_poly + pow_r(m) * dRb_vaddm
         enddo
c
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'a_H: ipar=',ipar
c-debug]
c
c__Skip over the parameters controlling the unused exponential cutoff:
c
      if ( ipar .gt. jpar ) then
         ipar = ipar + 1
         beta = parBMP(ipar)
         ipar = ipar + 1
         if ( beta .ne. 0.d0 ) then
            stop ' STOP -- Error: exp(beta*rho) term not implemented. '
c-noneed;            pow = parBMP(ipar)
c-noneed;            if ( pow .eq. 0.d0 ) then
c-noneed;               rhop = 1.d0
c-noneed;            else if ( pow .eq. 2.d0 ) then
c-noneed;               rhop = rhosq
c-noneed;            else
c-noneed;               rhop = rhosq**( 0.5d0 * pow )
c-noneed;            endif
c-noneed;            a_H_poly = a_H_poly * exp( min( - beta * rhop , hi_exp ) )
         endif
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'B_H: ipar=',ipar
c-debug]
c
c__Compute the F_H summations that depend on R_A or R_B alone (for V_H term):
c
      R_A_poly = 1.d0
      R_B_poly = 1.d0
c
      if ( ideriv .le. 0 ) then
c
         do j = -1, -JH, -1
            ipar = ipar + 1
            R_A_poly = R_A_poly + parBMP(ipar) * pow_Ra(j)
            R_B_poly = R_B_poly + parBMP(ipar) * pow_Rb(j)
         enddo
c
      else
c
         d_R_A_poly = 0.d0
         d_R_B_poly = 0.d0
c
         do j = -1, -JH, -1
            ipar = ipar + 1
            R_A_poly = R_A_poly + parBMP(ipar) * pow_Ra(j)
            R_B_poly = R_B_poly + parBMP(ipar) * pow_Rb(j)
            d_R_A_poly = d_R_A_poly + parBMP(ipar) * d_pow_Ra(j)
            d_R_B_poly = d_R_B_poly + parBMP(ipar) * d_pow_Rb(j)
         enddo
c
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'b_H: ipar=',ipar
c-debug]
c
c__Compute the V_H term and its derivatives:
c
c-noneed;      V_H = A_H * ( 1.d0
c-noneed;     $     - ( 1.d0 + a_H_poly ) * exp( min( - a * r , hi_exp ) ) )
c-noneed;     $     * ( exp( min( - dlambda_p * R_A , hi_exp ) ) * R_A_poly
c-noneed;     $     + exp( min( - dlambda_p * R_B , hi_exp ) ) * R_B_poly )
c
      f_H_Ra = exp( - dlambda_p * R_A ) * R_A_poly
      f_H_Rb = exp( - dlambda_p * R_B ) * R_B_poly
c
      a_H_fac = 1.d0 - ( 1.d0 + a_H_poly ) * exp( - a * r )
      V_H = A_H * a_H_fac * ( f_H_Ra + f_H_Rb )
c
      if ( ideriv .gt. 0 ) then
c
         dr_V_H = A_H * ( a * ( 1.d0 + a_H_poly ) - dr_a_H_poly )
     $        * exp( - a * r ) * ( f_H_Ra + f_H_Rb )
         dRa_V_H = A_H * ( a_H_fac * ( exp( - dlambda_p * R_A )
     $        * d_R_A_poly - dlambda_p * f_H_Ra )
     $        - dRa_A_H_poly * exp( - a * r ) * ( f_H_Ra + f_H_Rb ) )
         dRb_V_H = A_H * ( a_H_fac * ( exp( - dlambda_p * R_B )
     $        * d_R_B_poly - dlambda_p * f_H_Rb )
     $        - dRb_A_H_poly * exp( - a * r ) * ( f_H_Ra + f_H_Rb ) )
c
      endif
c
c-debug[
c-debug;      if ( idebug .gt. 0 .and. mod(idebug,modebug) .eq. 0 ) then
c-debug;         write(6,'(" ipar",i4,"  a_H_poly",f15.9,"  V_H",f15.9)')
c-debug;     $        ipar, a_H_poly, V_H
c-debug;      endif
c-debug]
c
c__Compute the summation for the V_cb term:
c
      A_cb_poly = 0.d0
c
      jpar = ipar
      ieq = 0
c
      if ( ideriv .le. 0 ) then
c
         do m = 0, Mrcb
            vaddm = 0.d0
            do L = 0, Lcb, 2
               vadd = 0.d0
               do i = Icb, Jcb
                  if ( i .ne. ieq ) then
                     ipar = ipar + 1
                     vadd = vadd + parBMP(ipar) * pow_Rbar(i)
                  endif
               enddo
               vaddm = vaddm + PL_eta(L) * vadd
               if ( L .eq. 2 ) ieq = 99
            enddo
            A_cb_poly = A_cb_poly + pow_r(m) * vaddm
            ieq = 99
         enddo
c
      else
c
         dr_A_cb_poly = 0.d0
         dRa_A_cb_poly = 0.d0
         dRb_A_cb_poly = 0.d0
c
         do m = 0, Mrcb
            vaddm = 0.d0
            dr_vaddm = 0.d0
            dRa_vaddm = 0.d0
            dRb_vaddm = 0.d0
            do L = 0, Lcb, 2
               vadd = 0.d0
               d_vadd = 0.d0
               do i = Icb, Jcb
                  if ( i .ne. ieq ) then
                     ipar = ipar + 1
                     vadd = vadd + parBMP(ipar) * pow_Rbar(i)
                     d_vadd = d_vadd + parBMP(ipar) * dRab_pow_Rbar(i)
                  endif
               enddo
               vaddm = vaddm + PL_eta(L) * vadd
               dr_vaddm = dr_vaddm + dr_PL_eta(L) * vadd
               dRa_vaddm = dRa_vaddm + dRa_PL_eta(L) * vadd
     $              + PL_eta(L) * d_vadd
               dRb_vaddm = dRb_vaddm + dRb_PL_eta(L) * vadd
     $              + PL_eta(L) * d_vadd
               if ( L .eq. 2 ) ieq = 99
            enddo
            A_cb_poly = A_cb_poly + pow_r(m) * vaddm
            dr_A_cb_poly = dr_A_cb_poly + d_pow_r(m) * vaddm
     $           + pow_r(m) * dr_vaddm
            dRa_A_cb_poly = dRa_A_cb_poly + pow_r(m) * dRa_vaddm
            dRb_A_cb_poly = dRb_A_cb_poly + pow_r(m) * dRb_vaddm
            ieq = 99
         enddo
c
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'acb: ipar=',ipar
c-debug]
c
c__Skip over the parameters controlling the unused exponential cutoff:
c
      if ( ipar .gt. jpar ) then
         ipar = ipar + 1
         beta = parBMP(ipar)
         ipar = ipar + 1
         if ( beta .ne. 0.d0 ) then
            stop ' STOP -- Error: exp(beta*rho) term not implemented. '
c-noneed;            pow = parBMP(ipar)
c-noneed;            if ( pow .eq. 0.d0 ) then
c-noneed;               rhop = 1.d0
c-noneed;            else if ( pow .eq. 2.d0 ) then
c-noneed;               rhop = rhosq
c-noneed;            else
c-noneed;               rhop = rhosq**( 0.5d0 * pow )
c-noneed;            endif
c-noneed;            A_cb_poly = A_cb_poly
c-noneed;     $           * exp( min( - beta * rhop , hi_exp ) )
         endif
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'Bcb: ipar=',ipar
c-debug]
c
c__Compute the Lambda_cb exponent, and the V_cb term and derivatives:
c
c-noneed;      dlambda_cb = dlcb_con
c-noneed;     $     + dlcb_norm * exp( min( dlcb_ex * r , hi_exp ) )
c-noneed;c
c-noneed;      V_cb = ( A_cb_poly + Acb_norm * ( 1.d0 + alpha_cb * P2_eta ) )
c-noneed;     $     * exp( min( Acb_ex * r - dlambda_cb * Rbar , hi_exp ) )
c
      cb_lam_ex_fac = dlcb_norm * exp( dlcb_ex * r )
      dlambda_cb = dlcb_con + cb_lam_ex_fac
c
      cb_V_ex_fac = exp( Acb_ex * r - dlambda_cb * Rbar )
      V_cb = ( A_cb_poly + Acb_norm * ( 1.d0 + alpha_cb * P2_eta ) )
     $     * cb_V_ex_fac
c
      if ( ideriv .gt. 0 ) then
c
         dr_V_cb = V_cb * ( Acb_ex
     $        - dlcb_ex * cb_lam_ex_fac * Rbar )
     $        + ( dr_A_cb_poly + Acb_norm * alpha_cb * dr_PL_eta(2) )
     $        * cb_V_ex_fac
         dRa_V_cb = -0.5d0 * dlambda_cb * V_cb
     $        + ( dRa_A_cb_poly + Acb_norm * alpha_cb * dRa_PL_eta(2) )
     $        * cb_V_ex_fac
         dRb_V_cb = -0.5d0 * dlambda_cb * V_cb
     $        + ( dRb_A_cb_poly + Acb_norm * alpha_cb * dRb_PL_eta(2) )
     $        * cb_V_ex_fac
c
         dr_V_P = dr_V_H + dr_V_cb
         dRa_V_P = dRa_V_H + dRa_V_cb
         dRb_V_P = dRb_V_H + dRb_V_cb
c
      endif
c
c-debug[
c-debug;      if ( idebug .gt. 0 .and. mod(idebug,modebug) .eq. 0 ) then
c-debug;         write(6,'(" ipar",i4,"  A_cb_poly",f15.9,"  V_cb",f15.9)')
c-debug;     $        ipar, A_cb_poly, V_cb
c-debug;      endif
c-debug]
c
c__By definition, V_P = V_H + V_cb:
c
      V_P = V_H + V_cb
c
c__Compute the long-range factors A_d and alpha_d (for the V_d term):
c
      alpha_d = ( f2_alpha_d * r + f1_alpha_d ) * r
     $     * fac_alpha_d * exp( fex_alpha_d * r )
      alpha_d_eta = alpha_d * P2_eta
      A_d = ( ( f2_A_d * r + f1_A_d ) * r + f0_A_d )
     $     * fac_A_d * exp( fex_A_d * r )
c
c__Compute the summation that will multiply 1/(Rbar**6+C_1**6) in V_d:
c
      A_d_poly = 1.d0 + alpha_d_eta
c
      if ( ideriv .le. 0 ) then
c
         do L = 4, Ld, 2
            if ( Mrd .gt. 0 ) then
               vadd = 0.d0
               do m = 1, Mrd
                  ipar = ipar + 1
                  vadd = vadd + parBMP(ipar) * pow_r(m)
               enddo
               ipar = ipar + 1
c-noneed;               A_d_poly = A_d_poly + vadd * PL_eta(L)
c-noneed;     $              * exp( min( parBMP(ipar) * r , hi_exp ) )
               A_d_poly = A_d_poly + vadd * PL_eta(L)
     $              * exp( parBMP(ipar) * r )
            endif
         enddo
c
      else
c
         dr_alpha_d = ( 2.d0 * f2_alpha_d * r + f1_alpha_d )
     $        * fac_alpha_d * exp( fex_alpha_d * r )
     $        + fex_alpha_d * alpha_d
         dr_A_d = ( 2.d0 * f2_A_d * r + f1_A_d )
     $        * fac_A_d * exp( fex_A_d * r )
     $        + fex_A_d * A_d
c
         dr_alpha_d_eta = dr_alpha_d * P2_eta + alpha_d * dr_PL_eta(2)
         dr_A_d_poly = dr_alpha_d_eta
         dRa_A_d_poly = alpha_d * dRa_PL_eta(2)
         dRb_A_d_poly = alpha_d * dRb_PL_eta(2)
c
         do L = 4, Ld, 2
            if ( Mrd .gt. 0 ) then
               vadd = 0.d0
               d_vadd = 0.d0
               do m = 1, Mrd
                  ipar = ipar + 1
                  vadd = vadd + parBMP(ipar) * pow_r(m)
                  d_vadd = d_vadd + parBMP(ipar) * d_pow_r(m)
               enddo
               ipar = ipar + 1
c-noneed;               A_d_poly = A_d_poly + vadd * PL_eta(L)
c-noneed;     $              * exp( min( parBMP(ipar) * r , hi_exp ) )
               vaddm = exp( parBMP(ipar) * r )
               A_d_poly = A_d_poly + vadd * PL_eta(L) * vaddm
               dr_A_d_poly = dr_A_d_poly
     $              + ( d_vadd * PL_eta(L) + vadd * dr_PL_eta(L)
     $              + vadd * PL_eta(L) * parBMP(ipar) ) * vaddm
               dRa_A_d_poly = dRa_A_d_poly
     $              + vadd * dRa_PL_eta(L) * vaddm
               dRb_A_d_poly = dRb_A_d_poly
     $              + vadd * dRb_PL_eta(L) * vaddm
            endif
         enddo
c		! compute derivatives BEFORE the modification of A_d_poly:
c
         dr_A_d_poly = dr_A_d_poly / ( pow_Rbar(Nmin) + C_1**Nmin )
     $        + C_3 * dr_alpha_d_eta / ( pow_Rbar(Ncont) + C_2**Ncont )
         d_vadd = A_d_poly * dRab_pow_Rbar(Nmin)
     $        / ( pow_Rbar(Nmin) + C_1**Nmin )**2
     $        + C_3 * ( 1.d0 + alpha_d_eta ) * dRab_pow_Rbar(Ncont)
     $        / ( pow_Rbar(Ncont) + C_2**Ncont )**2
         dRa_A_d_poly = dRa_A_d_poly / ( pow_Rbar(Nmin) + C_1**Nmin )
     $        + C_3 * alpha_d * dRa_PL_eta(2)
     $        / ( pow_Rbar(Ncont) + C_2**Ncont ) - d_vadd
         dRb_A_d_poly = dRb_A_d_poly / ( pow_Rbar(Nmin) + C_1**Nmin )
     $        + C_3 * alpha_d * dRb_PL_eta(2)
     $        / ( pow_Rbar(Ncont) + C_2**Ncont ) - d_vadd
c
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'ad0: ipar=',ipar
c-debug]
c
c__Get the 1/R**6 term and the first 1/R**8 term of V_d:
c
      A_d_poly = A_d_poly / ( pow_Rbar(Nmin) + C_1**Nmin )
     $     + C_3 * ( 1.d0 + alpha_d_eta )
     $     / ( pow_Rbar(Ncont) + C_2**Ncont )
c-debug[
c-debug;      if ( ideriv .gt. 0 .and. idebug .gt. 0 ) write(6,
c-debug;     $     '("&6. -A_d*A_d_poly",f25.15," dr",f25.15," +",f25.15,
c-debug;     $     " dRa,b",2f25.15)')
c-debug;     $     -A_d * A_d_poly, -A_d * dr_A_d_poly, -dr_A_d * A_d_poly,
c-debug;     $     -A_d * dRa_A_d_poly, -A_d * dRb_A_d_poly
c-debug]
c
c__Get the summation for any remaining parts of V_d:
c
      if ( ideriv .le. 0 ) then
c
         do n = Ncont, NRd, 2
            ilo = 0
            vl = 0.d0
            jpar = ipar
            do L = 0, Ld, 2
               if ( Mrd .ge. ilo ) then
                  vadd = 0.d0
                  do m = ilo, Mrd
                     ipar = ipar + 1
                     vadd = vadd + parBMP(ipar) * pow_r(m)
                  enddo
                  ipar = ipar + 1
c-noneed;                  vl = vl + vadd * PL_eta(L)
c-noneed;     $                 * exp( min( parBMP(ipar) * r , hi_exp ) )
                  vl = vl + vadd * PL_eta(L) * exp( parBMP(ipar) * r )
               endif
               ilo = 1
            enddo
            if ( ipar .gt. jpar ) then
               ipar = ipar + 1
               A_d_poly = A_d_poly
     $              + vl / ( pow_Rbar(n) + parBMP(ipar)**n )
            endif
         enddo
c
      else
c
         do n = Ncont, NRd, 2
            ilo = 0
            vl = 0.d0
            dr_vl = 0.d0
            dRa_vl = 0.d0
            dRb_vl = 0.d0
            jpar = ipar
            do L = 0, Ld, 2
               if ( Mrd .ge. ilo ) then
                  vadd = 0.d0
                  d_vadd = 0.d0
                  do m = ilo, Mrd
                     ipar = ipar + 1
                     vadd = vadd + parBMP(ipar) * pow_r(m)
                     d_vadd = d_vadd + parBMP(ipar) * d_pow_r(m)
                  enddo
                  ipar = ipar + 1
                  vaddm = exp( parBMP(ipar) * r )
                  vl = vl + vadd * PL_eta(L) * vaddm
                  dr_vl = dr_vl + ( d_vadd * PL_eta(L)
     $                 + vadd * ( dr_PL_eta(L)
     $                 + parBMP(ipar) * PL_eta(L) ) ) * vaddm
                  dRa_vl = dRa_vl + vadd * dRa_PL_eta(L) * vaddm
                  dRb_vl = dRb_vl + vadd * dRb_PL_eta(L) * vaddm
               endif
               ilo = 1
            enddo
            if ( ipar .gt. jpar ) then
               ipar = ipar + 1
               A_d_poly = A_d_poly
     $              + vl / ( pow_Rbar(n) + parBMP(ipar)**n )
               dr_A_d_poly = dr_A_d_poly
     $              + dr_vl / ( pow_Rbar(n) + parBMP(ipar)**n )
               dRa_A_d_poly = dRa_A_d_poly
     $              + dRa_vl / ( pow_Rbar(n) + parBMP(ipar)**n )
     $              - vl * dRab_pow_Rbar(n)
     $              / ( pow_Rbar(n) + parBMP(ipar)**n )**2
               dRb_A_d_poly = dRb_A_d_poly
     $              + dRb_vl / ( pow_Rbar(n) + parBMP(ipar)**n )
     $              - vl * dRab_pow_Rbar(n)
     $              / ( pow_Rbar(n) + parBMP(ipar)**n )**2
c-debug[
c-debug;               if ( ideriv .gt. 0 .and. idebug .gt. 0 ) write(6,
c-debug;     $              '(i2,". -A_d*vl/(Rbar^n+par**n)",f26.21,
c-debug;     $              " dr",f26.21," +",f26.21/
c-debug;     $              "      dRa",f26.21," +",f26.21,
c-debug;     $              " dRb",f26.21," +",f26.21)')
c-debug;     $              n, -A_d * vl / ( pow_Rbar(n) + parBMP(ipar)**n ),
c-debug;     $              -A_d * dr_vl / ( pow_Rbar(n) + parBMP(ipar)**n ),
c-debug;     $              -dr_A_d * vl / ( pow_Rbar(n) + parBMP(ipar)**n ),
c-debug;     $              -A_d * dRa_vl
c-debug;     $              / ( pow_Rbar(n) + parBMP(ipar)**n ),
c-debug;     $              A_d * vl * dRab_pow_Rbar(n)
c-debug;     $              / ( pow_Rbar(n) + parBMP(ipar)**n )**2,
c-debug;     $              -A_d * dRb_vl
c-debug;     $              / ( pow_Rbar(n) + parBMP(ipar)**n ),
c-debug;     $              A_d * vl * dRab_pow_Rbar(n)
c-debug;     $              / ( pow_Rbar(n) + parBMP(ipar)**n )**2
c-debug]
            endif
         enddo
c
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'a_d: ipar=',ipar
c-debug]
c
c__Get the V_d term and its derivatives, and the V_d2 term:
c
      V_d = - A_d * A_d_poly
c
c-debug[
c-debug;      if ( idebug .gt. 0 .and. mod(idebug,modebug) .eq. 0 ) then
c-debug;         write(6,'(" ipar",i4,"  A_d_poly",f15.9,"  V_d",f15.9)')
c-debug;     $        ipar, A_d_poly, V_d
c-debug;      endif
c-debug]
c
      if ( ideriv .le. 0 ) then
c
         vaddm = harkcal * A4m_geng / ( bohrangs / B4_geng )**6
         Vd_hhe1 = vaddm / ( pow_Ra(6) + c_heh**6 )
         Vd_hhe2 = vaddm / ( pow_Rb(6) + c_heh**6 )
         if ( NRheh .ge. Ncont ) then
            do n = Ncont, NRheh, 2
               jpar = ipar + 1
               ipar = ipar + 2
               Vd_hhe1 = Vd_hhe1
     $              + parBMP(jpar) / ( pow_Ra(n) + parBMP(ipar)**n )
               Vd_hhe2 = Vd_hhe2
     $              + parBMP(jpar) / ( pow_Rb(n) + parBMP(ipar)**n )
            enddo
         endif
c-noneed;         a_ex = c_ex * ( r**2 + d_ex**2 )
c-noneed;     $        * exp( min( e_ex * r , hi_exp ) )
         a_ex = c_ex * ( pow_r(2) + d_ex**2 ) * exp( e_ex * r )
         omf_ex = a_ex / sqrt( b_ex**2 + a_ex**2 )
         f_ex = 1.d0 - omf_ex
         V_d2 = f_ex * ( Vd_hhe1 + Vd_hhe2 )
c-debug[
c-debug;         if ( idebug .gt. 0 .and. mod(idebug,modebug) .eq. 0 ) then
c-debug;            write(6,'(" ipar",i4,"  Vd_hhe1,2",2f15.9,"  V_d",f15.9)')
c-debug;     $           ipar, Vd_hhe1, Vd_hhe2, V_d
c-debug;         endif
c-debug]
c
      else
c
         vaddm = harkcal * A4m_geng / ( bohrangs / B4_geng )**6
         Vd_hhe1 = vaddm / ( pow_Ra(6) + c_heh**6 )
         Vd_hhe2 = vaddm / ( pow_Rb(6) + c_heh**6 )
         d_Vd_hhe1 = - vaddm * d_pow_Ra(6)
     $        / ( pow_Ra(6) + c_heh**6 )**2
         d_Vd_hhe2 = - vaddm * d_pow_Rb(6)
     $        / ( pow_Rb(6) + c_heh**6 )**2
c-debug[
c-debug;         if ( ideriv .gt. 0 .and. idebug .gt. 0 ) write(6,
c-debug;     $        '(i2,":He.  VHe(Ra,b)",2f26.21,"  d_VHe",2f26.21)')
c-debug;     $        6, Vd_hhe1, Vd_hhe2, d_Vd_hhe1, d_Vd_hhe2
c-debug;         if ( ideriv .gt. 0 .and. idebug .gt. 0 ) write(6,
c-debug;     $        '(i2,":He.  VHe(Ra,b)",2f26.21,"  d_VHe",2f26.21)')
c-debug;     $        6, Vd_hhe1, Vd_hhe2,
c-debug;     $        -vaddm * 6 * R_A**5 / ( pow_Ra(6) + c_heh**6 )**2,
c-debug;     $        -vaddm * 6 * R_B**5 / ( pow_Rb(6) + c_heh**6 )**2
c-debug]
         if ( NRheh .ge. Ncont ) then
            do n = Ncont, NRheh, 2
               jpar = ipar + 1
               ipar = ipar + 2
               Vd_hhe1 = Vd_hhe1
     $              + parBMP(jpar) / ( pow_Ra(n) + parBMP(ipar)**n )
               Vd_hhe2 = Vd_hhe2
     $              + parBMP(jpar) / ( pow_Rb(n) + parBMP(ipar)**n )
               d_Vd_hhe1 = d_Vd_hhe1 - d_pow_Ra(n) * parBMP(jpar)
     $              / ( pow_Ra(n) + parBMP(ipar)**n )**2
               d_Vd_hhe2 = d_Vd_hhe2 - d_pow_Rb(n) * parBMP(jpar)
     $              / ( pow_Rb(n) + parBMP(ipar)**n )**2
c-debug[
c-debug;               if ( ideriv .gt. 0 .and. idebug .gt. 0 ) write(6,
c-debug;     $              '(i2,":He. term(Ra,b)",2f26.21," d_term",2f26.21)')
c-debug;     $              n,
c-debug;     $              parBMP(jpar) / ( pow_Ra(n) + parBMP(ipar)**n ),
c-debug;     $              parBMP(jpar) / ( pow_Rb(n) + parBMP(ipar)**n ),
c-debug;     $              -d_pow_Ra(n) * parBMP(jpar)
c-debug;     $              / ( pow_Ra(n) + parBMP(ipar)**n )**2,
c-debug;     $              -d_pow_Rb(n) * parBMP(jpar)
c-debug;     $              / ( pow_Rb(n) + parBMP(ipar)**n )**2
c-debug;               if ( ideriv .gt. 0 .and. idebug .gt. 0 ) write(6,
c-debug;     $              '(i2,":He. term(Ra,b)",2f26.21," d_term",2f26.21)')
c-debug;     $              n,
c-debug;     $              parBMP(jpar) / ( pow_Ra(n) + parBMP(ipar)**n ),
c-debug;     $              parBMP(jpar) / ( pow_Rb(n) + parBMP(ipar)**n ),
c-debug;     $              -parBMP(jpar) * n * R_A**(n-1)
c-debug;     $              / ( pow_Ra(n) + parBMP(ipar)**n )**2,
c-debug;     $              -parBMP(jpar) * n * R_B**(n-1)
c-debug;     $              / ( pow_Rb(n) + parBMP(ipar)**n )**2
c-debug]
            enddo
         endif
c-debug[
c-debug;         if ( ideriv .gt. 0 .and. idebug .gt. 0 ) write(6,
c-debug;     $        '("ALL:He. VHe(Ra,b)",2f26.21,"  d_VHe",2f26.21)')
c-debug;     $        Vd_hhe1, Vd_hhe2, d_Vd_hhe1, d_Vd_hhe2
c-debug]
         a_ex = c_ex * ( pow_r(2) + d_ex**2 ) * exp( e_ex * r )
         dr_a_ex = e_ex * a_ex + c_ex * d_pow_r(2) * exp( e_ex * r )
         vadd = b_ex**2 + a_ex**2
         vaddm = sqrt( vadd )
         omf_ex = a_ex / vaddm
         f_ex = 1.d0 - omf_ex
         dr_f_ex = ( a_ex**2 / vadd - 1.d0 ) * dr_a_ex / vaddm
c
         V_d2 = f_ex * ( Vd_hhe1 + Vd_hhe2 )
c
         dr_V_d = - dr_A_d * A_d_poly - A_d * dr_A_d_poly
         dRa_V_d = - A_d * dRa_A_d_poly
         dRb_V_d = - A_d * dRb_A_d_poly
c
         dr_V_d2 = dr_f_ex * ( Vd_hhe1 + Vd_hhe2 )
         dRa_V_d2 = f_ex * d_Vd_hhe1
         dRb_V_d2 = f_ex * d_Vd_hhe2 
c
c-debug[
c-debug;         if ( ideriv .gt. 0 .and. idebug .gt. 0 ) then
c-debug;            write(6,'(" a_ex=",f26.21," dr_a_ex=",f26.21,
c-debug;     $           " sqrt(a_ex^2+b_ex^2)=",f26.21," f_ex=",f26.21,
c-debug;     $           " dr_f_ex=",f26.21)')
c-debug;     $           a_ex, dr_a_ex, vaddm, f_ex, dr_f_ex
c-debug;            write(6,'(" V_d+V_D=",f26.21," =",f26.21,2(" +",f26.21))')
c-debug;     $           V_d + Vd2, V_d, f_ex * Vd_hhe1, f_ex * Vd_hhe2
c-debug;            write(6,'(" dr_V_Dd=",f26.21," =",f26.21,2(" +",f26.21))')
c-debug;     $           dr_V_d + dr_V_d2, dr_V_d2,
c-debug;     $           -dr_A_d * A_d_poly, -A_d * dr_A_d_poly
c-debug;            write(6,'("dRa_V_Dd=",f26.21," =",f26.21," +",f26.21)')
c-debug;     $           dRa_V_d + dRa_V_d2, dRa_V_d2, dRa_V_d
c-debug;            write(6,'("dRb_V_Dd=",f26.21," =",f26.21," +",f26.21)')
c-debug;     $           dRb_V_d + dRb_V_d2, dRb_V_d2, dRb_V_d
c-debug;         endif
c-debug]
      endif
c-debug[
c-debug;      if ( idebug .gt. 0 ) write(6,*) 'aHe: ipar=',ipar
c-debug]
c
c-noneed;      V_mbe = 0.d0
c
c__There should not be any MBE term:
c
      if ( Mrmbe .ge. 1 .and. Mhembe .ge. 1 .and. MBEord .ge. 2 ) then
c
         stop ' STOP -- Error: MBE term not implemented. '
c
c-noneed;         ipar = ipar + 1
c-noneed;         beta_mbe = parBMP(ipar)
c-noneed;         ipar = ipar + 1
c-noneed;         beta = parBMP(ipar)
c-noneed;         ipar = ipar + 1
c-noneed;         pow = parBMP(ipar)
c-noneed;c
c-noneed;         if ( beta_mbe .eq. 0.d0 ) then
c-noneed;            pr = r
c-noneed;            pa = R_A
c-noneed;            pb = R_B
c-noneed;         else
c-noneed;            pr = r * exp( min( - beta_mbe * r , hi_exp ) )
c-noneed;            pa = R_A * exp( min( - beta_mbe * R_A , hi_exp ) )
c-noneed;            pb = R_B * exp( min( - beta_mbe * R_B , hi_exp ) )
c-noneed;         endif
c-noneed;c
c-noneed;         pow_pr(0) = 1.d0
c-noneed;         pow_pa(0) = 1.d0
c-noneed;         pow_pb(0) = 1.d0
c-noneed;         pow_pr(1) = pr
c-noneed;         pow_pa(1) = pa
c-noneed;         pow_pb(1) = pb
c-noneed;c
c-noneed;         ihi = min( MBEord - 1 , Mrmbe )
c-noneed;         do i = 2, ihi
c-noneed;            pow_pr(i) = pow_pr(i-1) * pr
c-noneed;         enddo
c-noneed;c
c-noneed;         khi = min( MBEord - 1 , Mhembe )
c-noneed;         do k = 2, khi
c-noneed;            pow_pa(k) = pow_pa(k-1) * pa
c-noneed;            pow_pb(k) = pow_pb(k-1) * pb
c-noneed;         enddo
c-noneed;c
c-noneed;         do k = 1, khi
c-noneed;            do j = 0, min( k , MBEord - 1 - k )
c-noneed;               powjk = pow_pa(k) * pow_pb(j) + pow_pa(j) * pow_pb(k)
c-noneed;               do i = 1, min( Mrmbe , MBEord - j - k )
c-noneed;                  ipar = ipar + 1
c-noneed;                  V_mbe = V_mbe + parBMP(ipar) * pow_pr(i) * powjk
c-noneed;               enddo
c-noneed;            enddo
c-noneed;         enddo
c-noneed;c
c-noneed;         if ( beta .ne. 0.d0 ) then
c-noneed;            if ( pow .eq. 0.d0 ) then
c-noneed;               rhop = 1.d0
c-noneed;            else if ( pow .eq. 2.d0 ) then
c-noneed;               rhop = rhosq
c-noneed;            else
c-noneed;               rhop = rhosq**( 0.5d0 * pow )
c-noneed;            endif
c-noneed;            V_mbe = V_mbe * exp( min( - beta * rhop , hi_exp ) )
c-noneed;         endif
c
      endif
c
c__Combine all the terms to get Vint, and apply any required short-distance
c  extrapolation factors:
c
c-noneed;      Vint_0 = V_C + V_P + V_mbe + V_d + V_d2
c
      Vint_0 = V_C + V_P + V_d + V_d2
      Vadd = Vcl * fnmr + Vmr * omfnmr
      Vint = Vint_0 * fVint_Rab + Vadd * ofVint_Rab
      Vint_BMP = Vint
c
c__Get the total energy V:
c
      V = Vhh + Vint
c
c__Get the derivatives and return:
c
      if ( ideriv .gt. 0 ) then
c
         dr_Vint0 = dr_V_C + dr_V_P + dr_V_d + dr_V_d2
         dRa_Vint0 = dRa_V_C + dRa_V_P + dRa_V_d + dRa_V_d2
         dRb_Vint0 = dRb_V_C + dRb_V_P + dRb_V_d + dRb_V_d2
c
         dVhh_drhh = d_Vhh
         dVhhdr = d_Vhh
c
         dVint_drhh = Vint_0 * dfVint_drhh + Vadd * dofVint_drhh
     $        + ( dr_Vint0 * dr_drhh + dRa_Vint0 * dR_A_drhh
     $        + dRb_Vint0 * dR_B_drhh ) * fVint_Rab
     $        + ( ( Vcl - Vmr ) * dr_fnmr + dr_Vmr * omfnmr )
     $        * ofVint_Rab
         dVintdr = dVint_drhh
c
         dVdr = dVhh_drhh + dVint_drhh
         dVdRa = Vint_0 * dfVint_dRa + Vadd * dofVint_dRa
     $        + ( dr_Vint0 * dr_dRa + dRa_Vint0 * dR_A_dRa
     $        + dRb_Vint0 * dR_B_dRa ) * fVint_Rab
     $        + ( dR_A_Vcl * fnmr + dRa_Vmr * omfnmr ) * ofVint_Rab
         dVdRb = Vint_0 * dfVint_dRb + Vadd * dofVint_dRb
     $        + ( dr_Vint0 * dr_dRb + dRa_Vint0 * dR_A_dRb
     $        + dRb_Vint0 * dR_B_dRb ) * fVint_Rab
     $        + ( dR_B_Vcl * fnmr + dRb_Vmr * omfnmr ) * ofVint_Rab
c
      endif
c
c-debug[
c-debug;      if ( idebug .gt. 0 ) then
c-debug;         if ( mod(idebug,modebug) .eq. 0 ) then
c-debug;            write(6,'(" ipar",i4,"  NuseBMP",i4,":  Vint,V",2f15.9)')
c-debug;     $           ipar, NuseBMP, Vint, V
c-debug;            if ( ipar .ne. NuseBMP )
c-debug;     $           stop ' STOP -- Error: ipar .ne. NuseBMP. '
c-debug;         endif
c-debug;         idebug = idebug + 1
c-debug;         if ( idebug .gt. min(ndebug,24800) .and. ipar .ne. NuseBMP )
c-debug;     $        stop ' STOP -- Error: ipar .ne. NuseBMP. '
c-debug;         if ( idebug .gt. ndebug ) stop
c-debug;      else if ( ipar .ne. NuseBMP ) then
c-debug;         stop ' STOP -- Error: ipar .ne. NuseBMP. '
c-debug;      endif
c-debug]
c
      return
      end
c
c**********************************************************************
c
c  Values for the fitted BMP parameters (and a few other things):
c  ==============================================================
c
      block data heh2bmp_data
c     -----------------------
c
      implicit double precision (a-h,o-z)
c
      parameter ( r0_schw=3.5284882d0, DD_schw=0.160979391d0,
     $     C6_schw=6.499027d0, C8_schw=124.3991d0, C10_schw=3285.828d0,
     $     r02_schw=r0_schw*r0_schw, r04_schw=r02_schw*r02_schw,
     $     r06_schw=r04_schw*r02_schw , v6_schw=-9.163021333d0,
     $     v8_schw=642.104267d0, v10_schw=-60100.21337d0, nrschw=15 )
c
      common /heh2_bmp_schwenke_par/ a_schw(-1:nrschw),
     $     rsw_h2, rlo_h2, aa_h2, bb_h2, cc_h2
      save /heh2_bmp_schwenke_par/
c
      parameter ( Lmax=6, Mmax=4, Imax=4, Imin=-Imax, Nmin=6, Nmax=12,
     $     Ncont=Nmin+2, Iminm1=Imin-1, Jmax=5 )
c
      common /vheh2fit_poly_BMP/ Mrac, Lac, Iac, Jac,
     $     MrH, LH, IH, JH, Mrcb, Lcb, Icb, Jcb, NRd, Ld, Mrd, NRheh,
     $     Mrmbe, Mhembe, MBEord
      save /vheh2fit_poly_BMP/
c
      parameter ( nparBMP = 500, ndelBMP = 22 )
      parameter ( NcontBMP = ndelBMP + 1, NUMuseBMP = 113,
     $     NzBMP = NUMuseBMP + 1, nzeroBMP = nparBMP - NUMuseBMP )
c
      common /vheh2fit_par_BMP/ parBMP(nparBMP), NuseBMP
      save /vheh2fit_par_BMP/
      equivalence (parBMP(1),fac_A_d), (parBMP(2),fac_alpha_d),
     $     (parBMP(3),e_ex), (parBMP(4),d_ex),
     $     (parBMP(5),c_ex), (parBMP(6),b_ex),
     $     (parBMP(7),c_heh), (parBMP(8),C_1),
     $     (parBMP(9),C_2), (parBMP(10),C_3),
     $     (parBMP(11),Acb_norm), (parBMP(12),Acb_ex),
     $     (parBMP(13),dlcb_con), (parBMP(14),dlcb_norm),
     $     (parBMP(15),dlcb_ex), (parBMP(16),alpha_cb),
     $     (parBMP(17),A_H), (parBMP(18),a),
     $     (parBMP(19),dlambda_p), (parBMP(20),AC_norm),
     $     (parBMP(21),AC_ex), (parBMP(22),dlambda_c)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
c-debug[
c-debug;      common /vheh2fit_debug/ idebug, ndebug, modebug
c-debug;      save /vheh2fit_debug/
c-debug;      data idebug / 0 /, ndebug / 5 /, modebug / 1 /
c-debug]
c
c__Coefficients for the exp( polynomial ) in Schwenke's Vhh formula:
c
      data a_schw /
     $      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
c__Coefficients for the Vhh switchover to 1/r at small r:
c
      data rsw_h2 / 0.1193892052031558d0 /,
     $     rlo_h2 / 0.04515649888961240d0 /,
     $     aa_h2 / 0.42018214292507589d0 /,
     $     bb_h2 / 0.052365088838362515d0 /,
     $     cc_h2 / -0.0011823120379919859d0 /
c
c__Summation limits for the BMP surface formulae:
c
      data Mrac / 2 /, Lac / 2 /, Iac / 1 /, Jac / 5 /,
     $     MrH / 2 /, LH / 4 /, IH / 2 /, JH / 5 /,
     $     Mrcb / 2 /, Lcb / 4 /, Icb / -1 /, Jcb / 2 /,
     $     NRd / 8 /, Ld / 2 /, Mrd / 1 /, NRheh / 8 /,
     $     Mrmbe / -1 /, Mhembe / -1 /, MBEord / -1 /
c
c__Values for the 22 "basic" BMP surface parameters ("set 0"):
c
      data fac_A_d / 9.5516591661053818d-01 /,
     $     fac_alpha_d / 1.1490370307030446d+00 /
      data e_ex / -0.5891d0 /, d_ex / 6.d0 /,
     $     c_ex / 4.790521d0 /, b_ex / 6.d0 /
      data c_heh / 6.7876662396d0 /, C_1 / 4.3d0 /,
     $     C_2 / 5.05d0 /, C_3 / 3.8899009280953464d+01 /
      data Acb_norm / 1.0582959794485983d+02 /,
     $     Acb_ex / 7.5061006517216222d-01 /
      data dlcb_con / 3.1543676848688045d+00 /,
     $     dlcb_norm / 6.0488753691407231d+00 /,
     $     dlcb_ex / -1.0994065977015686d+00 /,
     $     alpha_cb / 5.9692607315813592d-02 /
      data A_H / 1.3246982020776207d+01 /,
     $     a / 7.0855446283386281d-01 /,
     $     dlambda_p / 2.0187445727733659d+00 /
      data AC_norm / 3.0421502781460488d+00 /,
     $     AC_ex / -2.9518270147913062d-01 /,
     $     dlambda_c / 1.5291696586348755d+00 /
c
c__Values for all the other BMP surface parameters:
c
      data (parBMP(i),i=NcontBMP,40) /
     $     -3.2245741278721733d+00, -3.3923915688808903d+00,
     $      3.0366190247398825d+00, -3.7458579571476691d-01,
     $      2.0058364264027850d-01,  4.4485330445898369d-01,
     $     -2.1515658796599302d-01,  0.0000000000000000d+00,
     $      2.0000000000000000d+00, -7.1379662898515903d+00,
     $      9.9197153023524614d+00, -3.4146406738211605d+00,
     $     -1.0267991316234168d+00,  5.8028434574410026d-01,
     $     -1.4720231016947809d+00,  1.3554644863456033d-01,
     $     -4.7585270109180605d-01,  1.4367505830547849d+00 /
c
      data (parBMP(i),i=41,60) /
     $     -4.7334520724766938d-01,  5.8062663873139175d-01,
     $     -1.6542155822260933d+00,  8.3051095943962638d-01,
     $      2.5588665689274999d+00, -8.7173209988343139d+00,
     $      5.4238304944886195d+00, -4.3957957653941560d-01,
     $      2.7744280038342228d-02, -1.3523686924206093d+00,
     $     -4.0081804384866515d-01,  9.6368237436319570d-01,
     $      2.3898717874166805d-01, -1.5450944359620435d+00,
     $      4.9731560422656305d+00, -2.8811449910371150d+00,
     $      1.4158786559150798d+00, -5.0876027733785829d+00,
     $      5.7490813892685084d+00, -3.2071418115997308d-01 /
c
      data (parBMP(i),i=61,80) /
     $      1.6575795612286364d+00, -1.6804142392307708d+00,
     $      0.0000000000000000d+00,  2.0000000000000000d+00,
     $     -2.7394834515659459d+00,  5.0804499843276725d+00,
     $     -5.4622258639825612d+00,  2.8935597281130856d+00,
     $     -5.8270493202587503d-01, -5.5078795269042899d+01,
     $      4.4320917570646650d+01, -3.1092730230812169d+01,
     $      5.2843594361500500d+01,  1.1600129588116692d+01,
     $     -1.9621689448693786d+01,  5.2900919827870531d+01,
     $     -1.0563209997064696d+02,  9.3284630315302877d+01,
     $     -4.0700642406249479d+00, -3.7052672514863647d+00 /
c
      data (parBMP(i),i=81,100) /
     $     -1.1346705968918137d+02,  5.1149810534447710d+01,
     $     -2.4871166169100647d+00, -9.3509311706637448d+01,
     $      1.6378977148427445d+01,  2.2086395299803186d+01,
     $     -1.6804990474796777d+00, -1.5853123070178475d+01,
     $      4.2638368107796731d+00, -3.5227430714860844d+01,
     $      5.9768890292206303d+00,  3.0364436387303375d+01,
     $     -1.3377093590413359d+01,  5.4999375817695695d-01,
     $      4.8222772834638758d-02,  3.8117938365375686d+01,
     $     -2.7665764465468698d+01,  4.8614292878569989d+00,
     $     -2.7649521198926752d-01,  5.2950814671088633d+00 /
c
      data (parBMP(i),i=101,NUMuseBMP) /
     $      6.5071927960341434d+00, -6.4462518717043993d-01,
     $     -1.9002875938543518d-01,  0.0000000000000000d+00,
     $      2.0000000000000000d+00,  2.0032220787100972d+01,
     $     -5.5038952296393710d+01, -1.0050049355200854d+00,
     $      1.5517530802208812d+01, -1.0000000000000000d+00,
     $      5.0000000000000000d+00, -5.2155365376656796d+01,
     $      5.0000000000000000d+00 /
c
      data (parBMP(i),i=NzBMP,nparBMP) / nzeroBMP * 0.d0 /
c
c__The total number of parameters used (actually, this includes a few
c  unused parameters and excludes some parameters treated as constants):
c
      data NuseBMP / NUMuseBMP /
c
c__Default values for error control and short-distance extrapolation:
c
      data Rabmin_err / 0.001d0 /, rmin_err / 0.001d0 /,
     $     Verror / 1.d36 /, level_err / 1 /, ifixclose / 0 /
c
      end
c
c**********************************************************************
c
c  The original Muchnick and Russek He-H2 surface (with Schwenke's Vhh):
c  =====================================================================
c
      subroutine heh2MR( rhh,Ra,Rb, ideriv, V, dVdr,dVdRa,dVdRb,
     $                   Vhh, dVhhdr, Vint, dVintdr )
c     ----------------------------------------------------------
c
c__heh2MR returns the HeH2 potential (with the Schwenke VH2 added)
c  of P.Muchnick and A.Russek (1994), J.Chem.Phys. 100, 4336-4346
c
      implicit double precision (a-h,o-z)
c
      call heh2MR_nD( rhh,Ra,Rb, ideriv,0, V, dVdr,dVdRa,dVdRb,
     $                Vhh, dVhhdr, Vint, dVintdr )
c
      return
      end
c
c**********************************************************************
c
c  The original Muchnick and Russek Vint:
c  ======================================
c
      subroutine VintMR(rhh,Ra,Rb,ideriv,Vint,dVintdr,dVdRa,dVdRb)
c     ------------------------------------------------------------
c
c__VintMR returns the HeH2 interaction potential (without VH2) of
c  P.Muchnick and A.Russek (1994), J.Chem.Phys. 100, 4336-4346
c
      implicit double precision (a-h,o-z)
c
      call VintMR_nD(rhh,Ra,Rb,ideriv,0,Vint,dVintdr,dVdRa,dVdRb)
c
      return
      end
c
c**********************************************************************
c
c  The "modMR" surface (from BMP):
c  ===============================
c
      subroutine heh2modMR( rhh,Ra,Rb, ideriv, V, dVdr,dVdRa,dVdRb,
     $                      Vhh, dVhhdr, Vint, dVintdr )
c     -------------------------------------------------------------
c
c__heh2modMR returns the modified Muchnick & Russek HeH2 interaction
c  potential (from BMP) with the Schwenke H2 potential added
c
      implicit double precision (a-h,o-z)
c
      call heh2MR_nD( rhh,Ra,Rb, ideriv,2, V, dVdr,dVdRa,dVdRb,
     $                Vhh, dVhhdr, Vint, dVintdr )
c
      return
      end
c
c**********************************************************************
c
c  Vint of the "modMR" surface:
c  ============================
c
      subroutine VintmodMR(rhh,Ra,Rb,ideriv,Vint,dVintdr,dVdRa,dVdRb)
c     ---------------------------------------------------------------
c
c__VintmodMR returns the modified Muchnick & Russek HeH2 interaction
c  potential (without VH2)
c
      implicit double precision (a-h,o-z)
c
      call VintMR_nD(rhh,Ra,Rb,ideriv,2,Vint,dVintdr,dVdRa,dVdRb)
c
      return
      end
c
c**********************************************************************
c
c  The MR surfaces as a function of the cartesian coordinates:
c  ===========================================================
c
      subroutine heh2MR_cc(ccH1,ccH2,ccHe,ideriv,nD,V,dVdr,dVdRa,dVdRb,
     $         Vhh,dVhhdr,Vint,dVintdr,r,Ra,Rb,dVdccH1,dVdccH2,dVdccHe)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
      dimension ccH1(3), ccH2(3), ccHe(3),
     $     dVdccH1(3), dVdccH2(3), dVdccHe(3)
c
      call VintMR_cc(ccH1,ccH2,ccHe,ideriv,nD,Vint,dVintdr,dVdRa,
     $               dVdRb,r,Ra,Rb,dVdccH1,dVdccH2,dVdccHe)
c
      if ( Vint .eq. Verror ) then
         dVdr = 0.d0
         dVhhdr = 0.d0
         V = Verror
         Vhh = Verror
         return
      endif
c
      call Vhh_heh2bmp(r,ideriv,Vhh,dVhhdr)
c
      V = Vint + Vhh
      dVdr = dVintdr + dVhhdr
c
      if ( ideriv .gt. 0 ) then
c
         dVdccH1(1) = dVdccH1(1) + dVhhdr * ( ccH1(1) - ccH2(1) ) / r
         dVdccH1(2) = dVdccH1(2) + dVhhdr * ( ccH1(2) - ccH2(2) ) / r
         dVdccH1(3) = dVdccH1(3) + dVhhdr * ( ccH1(3) - ccH2(3) ) / r
         dVdccH2(1) = dVdccH2(1) - dVhhdr * ( ccH1(1) - ccH2(1) ) / r
         dVdccH2(2) = dVdccH2(2) - dVhhdr * ( ccH1(2) - ccH2(2) ) / r
         dVdccH2(3) = dVdccH2(3) - dVhhdr * ( ccH1(3) - ccH2(3) ) / r
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  Vint of the MR surfaces as a function of the cartesian coordinates:
c  ===================================================================
c
      subroutine VintMR_cc(ccH1,ccH2,ccHe,ideriv,nD,Vint,dVintdr,dVdRa,
     $                     dVdRb,r,Ra,Rb,dVintdccH1,dVintdccH2,dVdccHe)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      dimension ccH1(3), ccH2(3), ccHe(3),
     $     dVintdccH1(3), dVintdccH2(3), dVdccHe(3)
c
      r = sqrt( ( ccH1(1) - ccH2(1) )**2
     $     + ( ccH1(2) - ccH2(2) )**2 + ( ccH1(3) - ccH2(3) )**2 )
      Ra = sqrt( ( ccH2(1) - ccHe(1) )**2
     $     + ( ccH2(2) - ccHe(2) )**2 + ( ccH2(3) - ccHe(3) )**2 )
      Rb = sqrt( ( ccH1(1) - ccHe(1) )**2
     $     + ( ccH1(2) - ccHe(2) )**2 + ( ccH1(3) - ccHe(3) )**2 )
c
      call VintMR_nD(r,Ra,Rb,ideriv,nD,Vint,dVintdr,dVdRa,dVdRb)
c
      if ( ideriv .gt. 0 ) then
c
         dVintdccH1(1) = dVintdr * ( ccH1(1) - ccH2(1) ) / r
     $        + dVdRb * ( ccH1(1) - ccHe(1) ) / Rb
         dVintdccH1(2) = dVintdr * ( ccH1(2) - ccH2(2) ) / r
     $        + dVdRb * ( ccH1(2) - ccHe(2) ) / Rb
         dVintdccH1(3) = dVintdr * ( ccH1(3) - ccH2(3) ) / r
     $        + dVdRb * ( ccH1(3) - ccHe(3) ) / Rb
         dVintdccH2(1) = - dVintdr * ( ccH1(1) - ccH2(1) ) / r
     $        + dVdRa * ( ccH2(1) - ccHe(1) ) / Ra
         dVintdccH2(2) = - dVintdr * ( ccH1(2) - ccH2(2) ) / r
     $        + dVdRa * ( ccH2(2) - ccHe(2) ) / Ra
         dVintdccH2(3) = - dVintdr * ( ccH1(3) - ccH2(3) ) / r
     $        + dVdRa * ( ccH2(3) - ccHe(3) ) / Ra
         dVdccHe(1) = - dVdRb * ( ccH1(1) - ccHe(1) ) / Rb
     $        - dVdRa * ( ccH2(1) - ccHe(1) ) / Ra
         dVdccHe(2) = - dVdRb * ( ccH1(2) - ccHe(2) ) / Rb
     $        - dVdRa * ( ccH2(2) - ccHe(2) ) / Ra
         dVdccHe(3) = - dVdRb * ( ccH1(3) - ccHe(3) ) / Rb
     $        - dVdRa * ( ccH2(3) - ccHe(3) ) / Ra
c
      else
c
         do i = 1, 3
            dVintdccH1(i) = 0.d0
            dVintdccH2(i) = 0.d0
            dVdccHe(i) = 0.d0
         enddo
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  The MR surfaces as a function of minimal cartesian coordinates:
c  ===============================================================
c
      subroutine heh2MR_cc3(a2,y3,z3,ideriv,nD, V,dVdr,dVdRa,dVdRb,Vhh,
     $                   dVhhdr,Vint,dVintdr,r,Ra,Rb,dVda2,dVdy3,dVdz3)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
      call VintMR_cc3(a2,y3,z3,ideriv,nD,Vint,dVintdr,dVdRa,
     $                dVdRb,r,Ra,Rb,dVintda2,dVdy3,dVdz3)
c
      if ( Vint .eq. Verror ) then
         dVdr = 0.d0
         dVhhdr = 0.d0
         V = Verror
         Vhh = Verror
         return
      endif
c
      call Vhh_heh2bmp(r,ideriv,Vhh,dVhhdr)
c
      V = Vint + Vhh
      dVdr = dVintdr + dVhhdr
      dVda2 = dVintda2 + dVhhdr * 2.d0
c
      return
      end
c
c**********************************************************************
c
c  Vint of the MR surfaces as a function of minimal cartesian coordinates:
c  =======================================================================
c
      subroutine VintMR_cc3(a2,y3,z3,ideriv,nD, Vint,dVintdr,dVdRa,
     $                      dVdRb,r,Ra,Rb,dVintda2,dVdy3,dVdz3)
c     -------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      r = 2.d0 * a2
      Ra = sqrt( ( z3 - a2 )**2 + y3**2 )
      Rb = sqrt( ( z3 + a2 )**2 + y3**2 )
c
      call VintMR_nD(r,Ra,Rb,ideriv,nD,Vint,dVintdr,dVdRa,dVdRb)
c
      if ( ideriv .gt. 0 ) then
c
         dr_da2 = 2.d0
         dRa_da2 = - ( z3 - a2 ) / Ra
         dRb_da2 = ( z3 + a2 ) / Rb
         dRa_dy3 = y3 / Ra
         dRb_dy3 = y3 / Rb
         dRa_dz3 = ( z3 - a2 ) / Ra
         dRb_dz3 = ( z3 + a2 ) / Rb
c
         dVintda2 = dVintdr * dr_da2
     $        + dVdRa * dRa_da2 + dVdRb * dRb_da2
         dVdy3 = dVdRa * dRa_dy3 + dVdRb * dRb_dy3
         dVdz3 = dVdRa * dRa_dz3 + dVdRb * dRb_dz3
c
      else
c
         dVintda2 = 0.d0
         dVdy3 = 0.d0
         dVdz3 = 0.d0
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  The MR surfaces as a function of angular coordinates:
c  =====================================================
c
      subroutine heh2MR_ang(rhh,rh2he,gammarad,ideriv,nD, V,dVdr,dVdRa,
     $      dVdRb,Vhh,dVhhdr,Vint,dVintdr,Ra,Rb,dVdrhh,dVdrh2he,dVdgam)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
      call VintMR_ang(rhh,rh2he,gammarad,ideriv,nD,Vint,dVintdr,
     $                dVdRa,dVdRb,Ra,Rb,dVintdrhh,dVdrh2he,dVdgam)
c
      if ( Vint .eq. Verror ) then
         dVdr = 0.d0
         dVhhdr = 0.d0
         V = Verror
         Vhh = Verror
         return
      endif
c
      call Vhh_heh2bmp(rhh,ideriv,Vhh,dVhhdr)
c
      V = Vint + Vhh
      dVdr = dVintdr + dVhhdr
      dVdrhh = dVintdrhh + dVhhdr
c
      return
      end
c
c**********************************************************************
c
c  Vint of the MR surfaces as a function of angular coordinates:
c  =============================================================
c
      subroutine VintMR_ang(rhh,rh2he,gammarad,ideriv,nD, Vint,dVintdr,
     $                     dVdRa,dVdRb,Ra,Rb,dVintdrhh,dVdrh2he,dVdgam)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      cgam = cos(gammarad)
      sgam = sin(gammarad)
      cgamR = cgam * rh2he
      sgamR = sgam * rh2he
      Ra = sqrt( ( cgamR - 0.5d0 * rhh )**2 + sgamR**2 )
      Rb = sqrt( ( cgamR + 0.5d0 * rhh )**2 + sgamR**2 )
c
      call VintMR_nD(rhh,Ra,Rb,ideriv,nD,Vint,dVintdr,dVdRa,dVdRb)
c
      if ( ideriv .gt. 0 ) then
c
         dRa_drhh = -0.5d0 * ( cgamR - 0.5d0 * rhh ) / Ra
         dRb_drhh = 0.5d0 * ( cgamR + 0.5d0 * rhh ) / Rb
         dRa_drh2he = ( cgam * ( cgamR - 0.5d0 * rhh )
     $        + sgam * sgamR ) / Ra
         dRb_drh2he = ( cgam * ( cgamR + 0.5d0 * rhh )
     $        + sgam * sgamR ) / Rb
         dRa_dgam = ( cgamR * sgamR
     $        - sgamR * ( cgamR - 0.5d0 * rhh ) ) / Ra
         dRb_dgam = ( cgamR * sgamR
     $        - sgamR * ( cgamR + 0.5d0 * rhh ) ) / Rb
c
         dVintdrhh = dVintdr + dVdRa * dRa_drhh + dVdRb * dRb_drhh
         dVdrh2he = dVdRa * dRa_drh2he + dVdRb * dRb_drh2he
         dVdgam = dVdRa * dRa_dgam + dVdRb * dRb_dgam
c
      else
c
         dVintdrhh = 0.d0
         dVdrh2he = 0.d0
         dVdgam = 0.d0
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  The MR surfaces as a function of the three distances:
c  =====================================================
c
      subroutine heh2MR_nD( rhh,Ra,Rb, ideriv,nD, V, dVdr,dVdRa,dVdRb,
     $                      Vhh, dVhhdr, Vint, dVintdr )
c     ----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
      call VintMR_nD(rhh,Ra,Rb,ideriv,nD,Vint,dVintdr,dVdRa,dVdRb)
c
      if ( Vint .eq. Verror ) then
         dVdr = 0.d0
         dVhhdr = 0.d0
         V = Verror
         Vhh = Verror
         return
      endif
c
      call Vhh_heh2bmp(rhh,ideriv,Vhh,dVhhdr)
c
      V = Vint + Vhh
      dVdr = dVintdr + dVhhdr
c
      return
      end
c
c**********************************************************************
c
c  Vint of the MR surfaces as a function of the three distances:
c  =============================================================
c
      subroutine VintMR_nD(rhh,Ra,Rb,ider,nD, Vint,dVintdr,dVdRa,dVdRb)
c     -----------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      parameter ( harkcal=1.d0/627.5096d0, bohrangs=0.52917725d0 )
      parameter ( A4m_geng=-3.597516D-2, B4_geng=3.2D0 )
c
      parameter ( dlambda_c = 4.3d0, dlambda_p = 1.940d0,
     $     A_H = 1.234d0, alpha_cb = 0.037d0, a = 1.935d0,
     $     C_1 = 4.3d0, C_2 = 5.05d0, C_3 = 30.45d0,
     $     C_1_six = C_1 * C_1 * C_1 * C_1 * C_1 * C_1,
     $     C_2_eight = C_2 * C_2 * C_2 * C_2 * C_2 * C_2 * C_2 * C_2 )
c
      parameter ( fac_alpha_d = 1.261d0, fac_A_d = 0.98234d0 )
c
      common /VintMR_parts_BMP/ V_C, V_P, V_d, V_d2, V_H, V_cb
      save /VintMR_parts_BMP/
c
      common /VintMR_deriv_BMP/
     $     dr_V_C, dRa_V_C, dRb_V_C, dr_V_P, dRa_V_P, dRb_V_P,
     $     dr_V_d, dRa_V_d, dRb_V_d, dr_V_d2, dRa_V_d2, dRb_V_d2,
     $     dr_V_H, dRa_V_H, dRb_V_H, dr_V_cb, dRa_V_cb, dRb_V_cb
      save /VintMR_deriv_BMP/
c
      dVintdr = 0.d0
      dVdRa = 0.d0
      dVdRb = 0.d0
c
      r = rhh
      R_A = Ra
      R_B = Rb
c
      call geomchk_heh2bmp( r, R_A, R_B, Vint, ierr )
c
      if ( ierr .ne. 0 ) return
c
      eta = ( R_A - R_B ) / r
      Rbar = ( R_A + R_B ) * 0.5d0
c
      etasq = eta**2
      P2_eta = 1.5d0 * etasq - 0.5d0
c
      if ( ider .gt. 0 ) then
         dr_etasq = -2.d0 * etasq / r
         dRa_etasq = 2.d0 * eta / r
         dRb_etasq = - dRa_etasq
         dr_P2_eta = 1.5d0 * dr_etasq
         dRa_P2_eta = 1.5d0 * dRa_etasq
         dRb_P2_eta = 1.5d0 * dRb_etasq
         dRab_Rbar = 0.5d0
      endif
c
      if ( nD .le. 0 ) then
c
         alpha_d = 0.112d0 + 0.06d0 * r
         alpha_d_eta = alpha_d * etasq
         A_d = 0.81d0 + 1.92d0 * r
c
         if ( ider .gt. 0 ) then
            dr_alpha_d = 0.06d0
            dr_A_d = 1.92d0
            dr_alpha_d_eta = dr_alpha_d * etasq + alpha_d * dr_etasq
            dRa_alpha_d_eta = alpha_d * dRa_etasq
            dRb_alpha_d_eta = alpha_d * dRb_etasq
         endif
c
      else
c
         alpha_d = ( 0.238248d0 * r - 0.126241d0 ) * r
     $        * fac_alpha_d * exp( -0.8075d0 * r )
         alpha_d_eta = alpha_d * P2_eta
         A_d = ( ( 4.790521d0 * r - 2.494444d0 ) * r + 3.008850d0 )
     $        * fac_A_d * exp( -0.5891d0 * r )
c
         if ( ider .gt. 0 ) then
            dr_alpha_d = ( 2.d0 * 0.238248d0 * r - 0.126241d0 )
     $           * fac_alpha_d * exp( -0.8075d0 * r )
     $           - 0.8075d0 * alpha_d
            dr_A_d = ( 2.d0 * 4.790521d0 * r - 2.494444d0 )
     $           * fac_A_d * exp( -0.5891d0 * r )
     $           - 0.5891d0 * A_d
            dr_alpha_d_eta = dr_alpha_d * P2_eta + alpha_d * dr_P2_eta
            dRa_alpha_d_eta = alpha_d * dRa_P2_eta
            dRb_alpha_d_eta = alpha_d * dRb_P2_eta
         endif
c
      endif
c
      V_d = - A_d * ( 1.d0 + alpha_d_eta )
     $     * ( 1.d0 / ( Rbar**6 + C_1_six )
     $     + C_3 / ( Rbar**8 + C_2_eight ) )
c
      if ( ider .gt. 0 ) then
         dr_V_d = - ( dr_A_d * ( 1.d0 + alpha_d_eta )
     $        + A_d * dr_alpha_d_eta )
     $        * ( 1.d0 / ( Rbar**6 + C_1_six )
     $        + C_3 / ( Rbar**8 + C_2_eight ) )
         dRa_V_d = A_d * ( 1.d0 + alpha_d_eta )
     $        * ( 6.d0 * dRab_Rbar * Rbar**5
     $        / ( Rbar**6 + C_1_six )**2
     $        + C_3 * 8.d0 * dRab_Rbar * Rbar**7
     $        / ( Rbar**8 + C_2_eight )**2 )
     $        - A_d * dRa_alpha_d_eta
     $        * ( 1.d0 / ( Rbar**6 + C_1_six )
     $        + C_3 / ( Rbar**8 + C_2_eight ) )
         dRb_V_d = A_d * ( 1.d0 + alpha_d_eta )
     $        * ( 6.d0 * dRab_Rbar * Rbar**5
     $        / ( Rbar**6 + C_1_six )**2
     $        + C_3 * 8.d0 * dRab_Rbar * Rbar**7
     $        / ( Rbar**8 + C_2_eight )**2 )
     $        - A_d * dRb_alpha_d_eta
     $        * ( 1.d0 / ( Rbar**6 + C_1_six )
     $        + C_3 / ( Rbar**8 + C_2_eight ) )
      endif
c
      if ( nD .le. 1 ) then
c
         V_d2 = 0.d0
         dr_V_d2 = 0.d0
         dRa_V_d2 = 0.d0
         dRb_V_d2 = 0.d0
c
      else
c
         Vd_hhe1 = harkcal * A4m_geng
     $        / ( ( bohrangs * R_A / 3.2d0 )**6 + 2.d0 )
         Vd_hhe2 = harkcal * A4m_geng
     $        / ( ( bohrangs * R_B / 3.2d0 )**6 + 2.d0 )
         a_ex = 4.790521d0 * ( r**2 + 36.d0 ) * exp( -0.5891d0 * r )
         f_ex = 1.d0 - a_ex / sqrt( 36.d0 + a_ex**2 )
         V_d2 = f_ex * ( Vd_hhe1 + Vd_hhe2 )
c
         if ( ider .gt. 0 ) then
            d_Vd_hhe1 = -6.d0 * harkcal * A4m_geng * bohrangs
     $           * ( bohrangs * R_A / 3.2d0 )**5 / ( 3.2d0
     $           * ( ( bohrangs * R_A / 3.2d0 )**6 + 2.d0 )**2 )
            d_Vd_hhe2 = -6.d0 * harkcal * A4m_geng * bohrangs
     $           * ( bohrangs * R_B / 3.2d0 )**5 / ( 3.2d0
     $           * ( ( bohrangs * R_B / 3.2d0 )**6 + 2.d0 )**2 )
            dr_a_ex = 2.d0 * 4.790521d0 * r * exp( -0.5891d0 * r )
     $           - 0.5891d0 * a_ex
            dr_f_ex = ( a_ex**2 / ( 36.d0 + a_ex**2 ) - 1.d0 )
     $           * dr_a_ex / sqrt( 36.d0 + a_ex**2 )
            dr_V_d2 = dr_f_ex * ( Vd_hhe1 + Vd_hhe2 )
            dRa_V_d2 = f_ex * d_Vd_hhe1
            dRb_V_d2 = f_ex * d_Vd_hhe2
         endif
c
      endif
c
      A_cb = 8.045d0 * exp( -0.360d0 * r )
c
      dlambda_cb = 1.388d0 + 1.150d0 * exp( -0.622d0 * r )
c
      V_cb = A_cb * ( 1.d0 + alpha_cb * P2_eta )
     $     * exp( - dlambda_cb * Rbar )
c
      V_H = A_H * ( 1.d0 - exp( - a * r ) )
     $     * ( exp( - dlambda_p * R_A ) + exp( - dlambda_p * R_B ) )
c
      V_P = V_H + V_cb
c
      A_C = 2.d0 - 4.34d0 * r * exp( -1.72d0 * r )
c
      V_C = A_C * ( 1.d0 + ( 2.d0 / A_C - 1.d0 ) * P2_eta )
     $     * ( exp( - dlambda_c * R_A ) / R_A
     $     + exp( - dlambda_c * R_B ) / R_B )
c
      Vint = V_C + V_P + V_d + V_d2
c
      if ( ider .gt. 0 ) then
c
         dr_A_cb = -0.360d0 * A_cb
         dr_dlambda_cb = -0.622d0 * 1.150d0 * exp( -0.622d0 * r )
c
         dr_V_cb = ( dr_A_cb * ( 1.d0 + alpha_cb * P2_eta )
     $        + A_cb * alpha_cb * dr_P2_eta )
     $        * exp( - dlambda_cb * Rbar )
     $        - dr_dlambda_cb * Rbar * V_cb
         dRa_V_cb = A_cb * alpha_cb * dRa_P2_eta
     $        * exp( - dlambda_cb * Rbar )
     $        - dlambda_cb * dRab_Rbar * V_cb
         dRb_V_cb = A_cb * alpha_cb * dRb_P2_eta
     $        * exp( - dlambda_cb * Rbar )
     $        - dlambda_cb * dRab_Rbar * V_cb
c
         dr_V_H = A_H * a * exp( - a * r )
     $        * ( exp( - dlambda_p * R_A ) + exp( - dlambda_p * R_B ) )
         dRa_V_H = - dlambda_p * A_H * ( 1.d0 - exp( - a * r ) )
     $        * exp( - dlambda_p * R_A )
         dRb_V_H = - dlambda_p * A_H * ( 1.d0 - exp( - a * r ) )
     $        * exp( - dlambda_p * R_B )
c
         dr_A_C = ( 1.72d0 * r - 1.d0 ) * 4.34d0 * exp( -1.72d0 * r )
c
         dr_V_C = ( dr_A_C * ( 1.d0 + ( 2.d0 / A_C - 1.d0 ) * P2_eta )
     $        + A_C * ( ( 2.d0 / A_C - 1.d0 ) * dr_P2_eta
     $        - 2.d0 * dr_A_C * P2_eta / A_C**2 ) )
     $        * ( exp( - dlambda_c * R_A ) / R_A
     $        + exp( - dlambda_c * R_B ) / R_B )
         dRa_V_C = A_C * ( ( 2.d0 / A_C - 1.d0 ) * dRa_P2_eta )
     $        * ( exp( - dlambda_c * R_A ) / R_A
     $        + exp( - dlambda_c * R_B ) / R_B )
     $        - ( dlambda_c + 1.d0 / R_A )
     $        * A_C * ( 1.d0 + ( 2.d0 / A_C - 1.d0 ) * P2_eta )
     $        * exp( - dlambda_c * R_A ) / R_A
         dRb_V_C = A_C * ( ( 2.d0 / A_C - 1.d0 ) * dRb_P2_eta )
     $        * ( exp( - dlambda_c * R_A ) / R_A
     $        + exp( - dlambda_c * R_B ) / R_B )
     $        - ( dlambda_c + 1.d0 / R_B )
     $        * A_C * ( 1.d0 + ( 2.d0 / A_C - 1.d0 ) * P2_eta )
     $        * exp( - dlambda_c * R_B ) / R_B
c
         dr_V_P = dr_V_H + dr_V_cb
         dRa_V_P = dRa_V_H + dRa_V_cb
         dRb_V_P = dRb_V_H + dRb_V_cb
c
         dVintdr = dr_V_d + dr_V_d2 + dr_V_P + dr_V_C
         dVdRa = dRa_V_d + dRa_V_d2 + dRa_V_P + dRa_V_C
         dVdRb = dRb_V_d + dRb_V_d2 + dRb_V_P + dRb_V_C
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  Schwenke's H-H interaction potential, with a fix-up at short distance:
c  ======================================================================
c
      subroutine Vhh_heh2bmp(r,ideriv,Vhh,dVhhdr)
c     -------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      parameter ( r0_schw=3.5284882d0, DD_schw=0.160979391d0,
     $     C6_schw=6.499027d0, C8_schw=124.3991d0, C10_schw=3285.828d0,
     $     r02_schw=r0_schw*r0_schw, r04_schw=r02_schw*r02_schw,
     $     r06_schw=r04_schw*r02_schw , v6_schw=-9.163021333d0,
     $     v8_schw=642.104267d0, v10_schw=-60100.21337d0, nrschw=15 )
c
      common /heh2_bmp_schwenke_par/ a_schw(-1:nrschw),
     $     rsw_h2, rlo_h2, aa_h2, bb_h2, cc_h2
      save /heh2_bmp_schwenke_par/
c
      dimension pow_r(nrschw)
c					! IF derivatives NOT required:
      if ( ideriv .le. 0 ) then
c
         dVhhdr = 0.d0
c					! usual r-case: use Schwenke formula:
         if ( r .gt. rsw_h2 ) then
c
            alphar_schw = a_schw(-1) / r + a_schw(0) + r * a_schw(1)
            pow_r(1) = r
c
            do m = 2, nrschw
               pow_r(m) = pow_r(m-1) * r
               alphar_schw = alphar_schw + pow_r(m) * a_schw(m)
            enddo
c
            ex_alph_schw = exp( alphar_schw )
            rr2 = pow_r(2) + r02_schw
            rr4 = pow_r(4) + r04_schw
            rr6 = pow_r(6) + r06_schw
            vsr_schw = DD_schw * ( ex_alph_schw - 1.d0 )**2 - DD_schw
            vlr_schw = - C6_schw / rr6 - C8_schw / rr4**2
     $           - C10_schw / rr2**5
c						! Vhh (from Schwenke formula)
            Vhh = vsr_schw + vlr_schw
c						! r < 0.1193892052031558 case:
         else
c						! Vhh: short-distance formulae
            if ( r .le. rlo_h2 ) then
               Vhh = 1.d0 / r
            else
               Vhh = ( ( cc_h2 / r + bb_h2 ) / r + aa_h2 ) / r
            endif
c
         endif
c					! ELSE: derivatives required:
      else
c					! usual r-case: use Schwenke formula:
         if ( r .gt. rsw_h2 ) then
c
            alphar_schw = a_schw(-1) / r + a_schw(0) + r * a_schw(1)
            d_alphar_schw = a_schw(1)
            pow_r(1) = r
c
            do m = 2, nrschw
               pow_r(m) = pow_r(m-1) * r
               alphar_schw = alphar_schw + pow_r(m) * a_schw(m)
               d_alphar_schw = d_alphar_schw
     $              + m * pow_r(m-1) * a_schw(m)
            enddo
c
            d_alphar_schw = d_alphar_schw - a_schw(-1) / pow_r(2)
c
            ex_alph_schw = exp( alphar_schw )
            rr2 = pow_r(2) + r02_schw
            rr4 = pow_r(4) + r04_schw
            rr6 = pow_r(6) + r06_schw
            vsr_schw = DD_schw * ( ex_alph_schw - 1.d0 )**2 - DD_schw
            vlr_schw = - C6_schw / rr6 - C8_schw / rr4**2
     $           - C10_schw / rr2**5
c						! Vhh (from Schwenke formula)
            Vhh = vsr_schw + vlr_schw
c								 ! d Vhh / d r
            d_vsr_schw = 2.d0 * DD_schw * ( ex_alph_schw - 1.d0 )
     $           * ex_alph_schw * d_alphar_schw
            d_vlr_schw = 6.d0 * C6_schw * pow_r(5) / rr6**2
     $           + 8.d0 * C8_schw * pow_r(3) / rr4**3
     $           + 10.d0 * C10_schw * r / rr2**6
            dVhhdr = d_vsr_schw + d_vlr_schw
c						! r < 0.1193892052031558 case:
         else
c						! Vhh: short-distance formulae
            if ( r .le. rlo_h2 ) then
               Vhh = 1.d0 / r
               dVhhdr = -1.d0 / r**2
            else
               Vhh = ( ( cc_h2 / r + bb_h2 ) / r + aa_h2 ) / r
               dVhhdr = ( ( -3.d0 * cc_h2 / r  - 2.d0 * bb_h2 ) / r
     $              - aa_h2 ) / r**2
            endif
c
         endif
c
      endif
c
      return
      end
c
c**********************************************************************
c
c  Check for invalid geometry or excessively short distance(s):
c  ============================================================
c
      subroutine geomchk_heh2bmp( r, R_A, R_B, V, ierr )
c     --------------------------------------------------
c
      implicit double precision (a-h,o-z)
c
      common /vheh2fit_fix_BMP/ Rabmin_err, rmin_err, Verror,
     $     level_err, ifixclose
      save /vheh2fit_fix_BMP/
c
      ierr = 0
      V = 0.d0
c
      if ( level_err .eq. 0 ) return
c
      if ( r .lt. rmin_err ) then
         ierr = 1
      else if ( R_A .lt. Rabmin_err ) then
         ierr = 2
      else if ( R_B .lt. Rabmin_err ) then
         ierr = 3
      else if ( max( r - R_A - R_B , R_A - r - R_B , R_B - r - R_A )
     $        .gt. 1.d-6 ) then
         ierr = 4
      endif
c
      if ( ierr .ne. 0 ) then
         V = Verror
         if ( level_err .lt. 0 ) stop
     $        ' STOP -- heh2bmp: Error in distances r, Ra, Rb. '
      endif
c
      return
      end
