SUBROUTINE HeHHp_DIM_pes(rs,E)
C SUBROUTINE TO CALCULATE THE DIM SURFACE FOR H2+ +He--->HeH+ +H
C  W.N. Whitton and P.J. Kurtz, JCP 64,3624(1976).
C    AB+C-->A+BC; A,B=H, C=He; R1=RH1H2, R2=RH2He, R3=RHeH1
C    dV/dR1, dV/dR2, and dV/dR3 are also calculated.
      IMPLICIT REAL*8(A-H,O-Z)
C     COMMON/POTCM/R,E,D
      DIMENSION RJ(3),R(3),DRJ(3),D(3),rs(1)
      DATA EAU/27.21161D0/
      DATA RAU/0.52917706/
      DATA U0/-2.8143D0/
      DATA IONE/1/,ITWO/2/
      DATA UNDFLO/1.D-60/
      U123=0.D0
      r(1)=rs(1)
      r(2)=rs(2)
      r(3)=rs(3)
      DO 100 J=1,3
      JJ=J
      R(J)=R(J)*RAU
      CALL HeHHp_DIM_V12(IONE,JJ,R(J),V1,DV1)
      CALL HeHHp_DIM_V12(ITWO,JJ,R(J),V2,DV2)
      Q=0.5D0*(V1+V2)
      D(J)=0.5D0*(DV1+DV2)
      RJ(J)=0.5D0*(V1-V2)
      DRJ(J)=0.5D0*(DV1-DV2)
      U123=U123+Q
      R(J)=R(J)/RAU
  100 CONTINUE
      RJSQ=RJ(1)**2+(RJ(2)-RJ(3))**2
      RJSQ=DSQRT(RJSQ)
      IF(RJSQ.LT.UNDFLO) RJSQ=UNDFLO
      RJSQI=1.D0/RJSQ
      U123=U123-RJSQ
      U123=U123-U0
      D(1)=D(1)-RJSQI*RJ(1)*DRJ(1)
      D(1)=D(1)*RAU/EAU
      D(2)=D(2)-RJSQI*(RJ(2)-RJ(3))*(DRJ(2))
      D(2)=D(2)*RAU/EAU
      D(3)=D(3)+RJSQI*(RJ(2)-RJ(3))*(DRJ(3))
      D(3)=D(3)*RAU/EAU
      U123=U123/EAU
      E=U123
c     D(1)=D(1)-1.D0/(R(1)**2)
c     D(2)=D(2)-1.D0/(R(2)**2)
c     D(3)=D(3)-1.D0/(R(3)**2)
      RETURN
      END
      SUBROUTINE HeHHp_DIM_V12(I,J,R,V,DVJ)
C Calculates V1 and V2 for the DIM method. Derivatives are also
C   calculated.
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION D(2,3),R0(2,3),BETA(2,3),GAMMA(2,3),RLAMDA(2,3),P(2)
      DATA D/2.8141D0,2.9014D0,2.0389D0,1.6431D0,2.0389D0,1.6431D0/
      DATA R0/1.0490D0,1.0680D0,0.7761D0,0.7226D0,0.7761D0,0.7226D0/
      DATA BETA/1.4049D0,1.4968D0,2.7936D0,3.0664D0,2.7936D0,3.0664D0/
      DATA GAMMA/-0.1374D0,-0.0573D0,0.0656D0,-0.0227D0,0.0656D0,
     ?-0.0227D0/
      DATA RLAMDA/0.0514D0,0.0015D0,0.0453D0,0.0681D0,0.0453D0,0.0681D0/
      DATA P/-1.D0,1.D0/
      RD=R-R0(I,J)
      B=BETA(I,J)*(1.D0+GAMMA(I,J)*RD+RLAMDA(I,J)*RD**2)
      DBJ=BETA(I,J)*(GAMMA(I,J)+2.D0*RLAMDA(I,J)*RD)
      VX=P(I)*DEXP(-B*RD)
      VX1=1.D0+VX
      V=D(I,J)*VX1**2-D(I,J)
      DVJ=-2.D0*D(I,J)*VX1*VX*(B+RD*DBJ)
      RETURN
      END