cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc
      subroutine HHH_DMBE_init()
      call prepot_HHH_DMBE()
      return
      end
cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc
      subroutine HHH_DMBE_pes(rs, v)
      implicit none
      real*8 rs(3), v
      !! Common block
      real*8 X(3),E,E1,E2,E3
      COMMON /POTCM_HHH_DMBE/ X,E,E1,E2,E3

      x(1)=rs(1)
      x(2)=rs(2)
      x(3)=rs(3)
      call pot_HHH_DMBE
      v=E

      return
      end
cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc

C   User must call prepot once for initialization of constants
C   Then call pot when you want the energy
C   Pass parameters through potcm
      SUBROUTINE PREPOT_HHH_DMBE
C***********************************************************************
C     Calculates double many-body expansion of the H3 potential.
C        POTE1 - energy of surface 1.
C        DH3L  - derivatives for surface 1.
C        POTE2 - energy of surface 2.
C        DH3U  - derivatives for surface 2.
C***********************************************************************
C          A.J.C. Varandas, F.B. Brown, C.A.Mead,D.G.Truhlar,and N.C.
C              Blais, J. Chem. Phys. 86, 6258 (1987).
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /POTCM_HHH_DMBE/ R(3), POTE1, DH3L(3)
      COMMON /POT2CM_HHH_DMBE/ POTE2, DH3U(3), COUP(3)
C  Common block for H3 potential parameters, set in BLOCK DATA H3
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
C  Common block for coordinates (computed in H3COOR)
      COMMON /H3COCM_HHH_DMBE/ PER, PER2, R12, R22, R32, QCOORD, DQ1, 
     $DQ2, DQ3,
     *   RHO, DRHO1, DRHO2, DRHO3, S, S2, DS1, DS2, DS3, CPHI3, DCPHI1,
     *   DCPHI2, DCPHI3
C  Common block for 2-body correlation energies and damping factors
C     (computed in H3COR2)
      COMMON /H3CRCM_HHH_DMBE/ CORRS(3), DCORRS(3), CORRT(3), DCORRT(3),
     *   DAMP(3,3), DDAMP(3,3)
      DIMENSION DLEP(3), DLEP2(3), DV2(3), DV3(3), DVA(3), DC3(3)
C
      DASY = 0.174474112D0
      NEXP = 4
      DO 10 ID = 1,3
         NEXP = NEXP + 2
         DAMPA(ID) = ALPH0/FLOAT(NEXP)**ALPH1
         DAMPB(ID) = BET0*EXP(-BET1*FLOAT(NEXP))
   10 CONTINUE
      RETURN
C
      ENTRY POT_HHH_DMBE
C
C Set up symmetry coordinates and their derivatives
      CALL H3COOR_HHH_DMBE (R(1), R(2), R(3))
C
C Set up 2-body correlation energies, dispersion damping terms, and
C    their derivatives
      CALL H3COR2_HHH_DMBE (R)
C
C Get Leps potentials and derivatives
      CALL H3LEPS_HHH_DMBE (R, VLEPS, VLEPS2, DLEP, DLEP2)
C
C Get VA term and its derivatives
      CALL H3VA_HHH_DMBE (R(1), R(2), R(3), VA, DVA)
C
C Get VII term and its derivatives
      CALL H3VII_HHH_DMBE (R(1), R(2), R(3), VII, DV2)
C
C Get VIII term and its derivatives
      CALL H3VIII_HHH_DMBE (R(1), R(2), R(3), VIII, DV3)
C
C Get 3-body correlation and its derivative
      CALL H3COR3_HHH_DMBE (R, CE3, DC3)
C
C 2-body correlation term
      CE2 = CORRS(1) + CORRS(2) + CORRS(3)
C
      VC = VA + VII + VIII + CE2 + CE3 + DASY
      POTE1 = VLEPS + VC
      POTE2 = VLEPS2 + VC
      DO 20 I = 1,3
         T = DVA(I) + DV2(I) + DV3(I) + DCORRS(I) + DC3(I)
         DH3L(I) = DLEP(I) + T
         DH3U(I) = DLEP2(I) + T
   20 CONTINUE
      RETURN
      END
      BLOCK DATA H3_HHH_DMBE
C***********************************************************************
C     Data for the DMBE H3 surface.
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
C
      DATA RHOL /2.4848D0/
      DATA ALPH0, ALPH1, BET0, BET1 /25.9528D0, 1.1868D0, 15.7381D0,
     *   0.09729D0/
      DATA ALPHA5, BETA1, BETA2, BETA3 / 8.2433033D-03, 0.53302897D0,
     *   0.39156612D-01, 0.69996945D0/
      DATA ALPH2 /4.735364D-1/
      DATA AL0, AL1, AL2, AL3 / 0.45024472D+01, -0.62467617D+01,
     *   0.40966542D+01, 0.21813012D+01/
      DATA AZ2 /4.453649D-4/
      DATA CD0, CD1 /6.333404D-03, -1.726839D-03/
      DATA CHH /  6.499027D0,      1.243991D+02,  3285.8D0/
      DATA CK0 / -1.372843D-01, -1.638459D-01, -1.973814D-01/
      DATA CK1 /  1.011204D0,      9.988099D-01,  9.399411D-01/
      DATA HFD, HFA1, HFA2, HFA3, HFGAM /0.15796326D0, 2.1977034D0,
     *   1.2932502D0, 0.64375666D0, 2.1835071D0/
      DATA H2RM, H2R0, H2RMT /1.401D0, 6.928203D0, 7.82D0/
      DATA H2TA, H2TB1, H2TB2, H2TB3, H2TB4 /0.448467D0, -0.056687D0,
     *   0.246831D0, -0.018419D0, 0.000598D0/
      DATA PI /3.14159265358979D0/
      DATA SQRT3 /1.73205080756887D0/
      DATA XPAR /-0.9286579D-02,  0.2811592D-03, -0.4665659D-05,
     *            0.2069800D-07,  0.2903613D+02, -0.2934824D+01,
     *            0.7181886D0,     -0.3753218D0,     -0.1114538D+01,
     *            0.2134221D+01, -0.4164343D0,      0.2022584D0,
     *           -0.4662687D-01, -0.4818623D+02,  0.2988468D0/
      END
      SUBROUTINE H3COOR_HHH_DMBE (R1, R2, R3)
C**********************************************************************
C   Calculates D3H symmetry coordinates and derivatives
C**********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
      COMMON /H3COCM_HHH_DMBE/ PER, PER2, R12, R22, R32, QCOORD, DQ1, 
     $DQ2, DQ3,
     *   RHO, DRHO1, DRHO2, DRHO3, S, S2, DS1, DS2, DS3, CPHI3, DCPHI1,
     *   DCPHI2, DCPHI3
C
      PER = R1+R2+R3
      PER2 = PER*PER
C   QCOORD and its derivatives
      R12 = R1*R1
      R22 = R2*R2
      R32 = R3*R3
      QCOORD = R12+R22+R32
      DQ1 = 2.0D0*R1
      DQ2 = 2.0D0*R2
      DQ3 = 2.0D0*R3
C   RHO and its derivatives
      RHO = SQRT(QCOORD/3.0D0)
      T = 1.0D0/(6.0D0*RHO)
      DRHO1 = T*DQ1
      DRHO2 = T*DQ2
      DRHO3 = T*DQ3
C   S, CPHI3 (cos(phi3)), and their derivatives
      GAMMA = 2.0D0*R12-R22-R32
      GAM2 = GAMMA*GAMMA
      DGM1 = 2.0D0*DQ1
      DGM2 = -DQ2
      DGM3 = -DQ3
      BETA = SQRT3*(R22-R32)
      BET2 = BETA*BETA
      DBT1 = 0.0D0
      DBT2 = SQRT3*DQ2
      DBT3 = -SQRT3*DQ3
      T12 = BET2 + GAM2
      T1 = SQRT(T12)
      S = T1/QCOORD
      S2 = S*S
      IF (S .EQ. 0.0D0) THEN
         DS1 = 0.0D0
         DS2 = 0.0D0
         DS3 = 0.0D0
C   For S=0, CPHI3 and its derivative should not be used anywhere but
C      set to zero anyway.
         CPHI3 = 0.0D0
         DCPHI1 = 0.0D0
         DCPHI2 = 0.0D0
         DCPHI3 = 0.0D0
      ELSE
         DS1 = S*((BETA*DBT1+GAMMA*DGM1)/T12 - DQ1/QCOORD)
         DS2 = S*((BETA*DBT2+GAMMA*DGM2)/T12 - DQ2/QCOORD)
         DS3 = S*((BETA*DBT3+GAMMA*DGM3)/T12 - DQ3/QCOORD)
         T2 = 1.0D0/(T1*T12)
         CPHI3 = GAMMA*(3.0D0*BET2 - GAM2)*T2
         T3 = 3.0D0*BETA*(3.0D0*GAM2 - BET2)*T2/T12
         DCPHI1 = T3*(GAMMA*DBT1 - BETA*DGM1)
         DCPHI2 = T3*(GAMMA*DBT2 - BETA*DGM2)
         DCPHI3 = T3*(GAMMA*DBT3 - BETA*DGM3)
      END IF
      RETURN
      END
      SUBROUTINE H3COR2_HHH_DMBE (RR)
C***********************************************************************
C     Calculates 2-body correlation energies for singlet and triplet
C     states of H2, the damping factors for the dispersion terms, and
C     their 1st derivatives
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION RR(3)
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
      COMMON /H3CRCM_HHH_DMBE/ CORRS(3), DCORRS(3), CORRT(3), DCORRT(3),
     *   DAMP(3,3), DDAMP(3,3)
C
C   Loop over three coordinates
      DO 10 J = 1,3
         R = RR(J)
         CORRS(J) = 0.0D0
         DCORRS(J) = 0.0D0
         CORRT(J) = 0.0D0
         DCORRT(J) = 0.0D0
         T1 = 1.0D0/R
         T = T1*T1
         T2 = T*T
         T1 = T2*T1
         NEXP = 4
C      Loop over terms in dispersion expansion
         DO 1 ID = 1,3
            NEXP = NEXP + 2
            FNEXP = FLOAT(NEXP)
            T2 = T2*T
            T1 = T1*T
C     singlet
C        Calculate damping function for the nth dispersion coefficient
C        and its 1st derivative
            DENOM = H2RM + 2.5D0*H2R0
            X = 2.0D0*R/DENOM
            TT = DAMPB(ID)*X
            TT1 = DAMPA(ID) + TT
            POL = X*TT1
            DPOL = TT1 + TT
            TT = EXP(-POL)
            TT1 = 1.0D0 - TT
            TT2 = TT1**(NEXP-1)
            D = TT1*TT2
            DD = 2.0D0*FNEXP*DPOL*TT*TT2/DENOM
C     store damping factors (including CHH coeff and 1/R**NEXP) for
C        later use in computing the 3-body correlation energy.
            DAMP(ID,J) = CHH(ID)*D*T2
            DDAMP(ID,J) = CHH(ID)*(DD*T2 - FNEXP*D*T1)
            CORRS(J) = CORRS(J) - DAMP(ID,J)
            DCORRS(J) = DCORRS(J) - DDAMP(ID,J)
C     triplet
C        Calculate damping function for the nth dispersion coefficient
C        and its 1st derivative
            DENOM = H2RMT + 2.5D0*H2R0
            X = 2.0D0*R/DENOM
            TT = DAMPB(ID)*X
            TT1 = DAMPA(ID) + TT
            POL = X*TT1
            DPOL = TT1 + TT
            TT = EXP(-POL)
            TT1 = 1.0D0 - TT
            TT2 = TT1**(NEXP-1)
            D = TT1*TT2
            DD = 2.0D0*FNEXP*DPOL*TT*TT2/DENOM
            CORRT(J) = CORRT(J) - CHH(ID)*D*T2
            DCORRT(J) = DCORRT(J) - CHH(ID)*(DD*T2 - FNEXP*D*T1)
    1    CONTINUE
   10 CONTINUE
      RETURN
      END
      SUBROUTINE H3COR3_HHH_DMBE (R, CE3, DC3)
C***********************************************************************
C     Calculates 3-body correlation energy and its 1st derivatives
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(3), DC3(3)
      DIMENSION G(3), GD(3), H(3), HD(3)
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
      COMMON /H3CRCM_HHH_DMBE/ CORRS(3), DCORRS(3), CORRT(3), DCORRT(3),
     *   DAMP(3,3), DDAMP(3,3)
C
      CE3 = 0.0D0
      DC3(1) = 0.0D0
      DC3(2) = 0.0D0
      DC3(3) = 0.0D0
C   Loop over terms in dispersion expansion; dispersion damping factors
C      are computed in H3COR2 and passed through COMMON /H3CRCM/.
      DO 30 ID = 1,3
C   Loop over 3 coordinates
         DO 10 I = 1,3
            RR = R(I)
C     Calculate G function and its 1st derivative
            T = CK0(ID)*EXP(-CK1(ID)*(RR-H2RM))
            G(I) = 1.0D0 + T
            GD(I) = -CK1(ID)*T
C     Calculate H function and its 1st derivative.
            T = CK1(ID)*RR
            SGNT = 1.0D0
            IF (T .LT. 0.0D0) THEN
               T = -T
               SGNT = -1.0D0
            END IF
            T = EXP(-T)
            T2 = T*T
            T1 = 1.0D0/(1.0D0+T2)
            HYSEC = 2.0D0*T*T1
            HYTAN = SGNT*(1.0D0-T2)*T1
            T1 = HYTAN**5
            H(I) = HYTAN*T1
            HD(I) = 6.0D0*CK1(ID)*T1*HYSEC*HYSEC
   10    CONTINUE
         DO 20 I = 1,3
            I2 = MOD(I,3) + 1
            I3 = MOD(I+1,3) + 1
            T = 1.0D0 - 0.5D0*(G(I2)*H(I3) + G(I3)*H(I2))
            CE3 = CE3 + T*DAMP(ID,I)
            DC3(I) = DC3(I) + T*DDAMP(ID,I)
            DC3(I2) = DC3(I2) -
     *         0.5D0*(GD(I2)*H(I3)+G(I3)*HD(I2))*DAMP(ID,I)
            DC3(I3) = DC3(I3) -
     *         0.5D0*(G(I2)*HD(I3)+GD(I3)*H(I2))*DAMP(ID,I)
   20    CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE H3LEPS_HHH_DMBE (R, VLEPS, VLEPS2, DLEP, DLEP2)
C***********************************************************************
C     Calculates 3-body extended-Hartree-Fock energy defined by a
C     LEPS-type function.
C***********************************************************************
C   VLEPS  IS LEPS LOWER SURFACE.
C   VLEPS2 IS LEPS UPPER SURFACE.
C   DLEP(3) ARE LEPS LOWER SURFACE DERIVATIVES
C   DLEP2(3) "   "   UPPER    "       "
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(3), DLEP(3), DLEP2(3)
      DIMENSION DF(3), DQ(3), XJ(3), DXJ(3,3)
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
      COMMON /H3CRCM_HHH_DMBE/ CORRS(3), DCORRS(3), CORRT(3), DCORRT(3),
     *   DAMP(3,3), DDAMP(3,3)
      COMMON /H3COCM_HHH_DMBE/ PER, PER2, R12, R22, R32, QCOORD, DQ1, 
     $DQ2, DQ3,
     *   RHO, DRHO1, DRHO2, DRHO3, S, S2, DS1, DS2, DS3, CPHI3, DCPHI1,
     *   DCPHI2, DCPHI3
C
C  Calculate switching term for the Hartree-Fock component of the
C     diatomic triplet state function.  The notation of Thompson et al.
C     (J. Chem. Phys., 82,5597,1985; and references therein) is used
C     throughout.
C
      T1 = -AZ2*QCOORD
      IF (S .EQ. 0.0D0) THEN
         F = EXP(T1*QCOORD)
         T1 = 2.0D0*T1*F
         DF(1) = T1*DQ1
         DF(2) = T1*DQ2
         DF(3) = T1*DQ3
      ELSE
         T2 = 1.0D0 + S*S2*CPHI3
         F = EXP(T1*QCOORD*T2)
         T1 = T1*F
         TDQ = 2.0D0*T2
         T2 = S2*QCOORD
         TDS = 3.0D0*CPHI3*T2
         TDC = S*T2
         DF(1) = T1*(TDQ*DQ1 + TDS*DS1 + TDC*DCPHI1)
         DF(2) = T1*(TDQ*DQ2 + TDS*DS2 + TDC*DCPHI2)
         DF(3) = T1*(TDQ*DQ3 + TDS*DS3 + TDC*DCPHI3)
      END IF
      OMF = 1.0D0 - F
C
      QSUM = 0.0D0
      DO 10 I = 1,3
         DQ(I) = 0.0D0
   10 CONTINUE
C  Loop over 3 coordinates
      DO 50 I = 1,3
         RR = R(I)
         RRI = 1.0D0/RR
C     Calculate extended-Hartree-Fock curve for ground singlet state of
C     H2 [A.J.C. Varandas & J.D. Silva, J. Chem. Soc. Faraday II
C     (submitted)] and 1st derivative of 2-body extended-Hartree-Fock
C     curve for the ground-singlet state of H2.
         DR = RR - H2RM
         T = EXP(-HFGAM*DR)
         ES = -HFD*(1.0D0 + DR*(HFA1 + DR*(HFA2 + DR*HFA3)))*T
         DESDR = -HFGAM*ES - HFD*(HFA1 + DR*(2.0D0*HFA2 +
     *      3.0D0*DR*HFA3))*T
C     Compute the HFACE (Hartree-Fock-approximate correlation energy)
C     potential for the lowest triplet state of H2 [A.J.C. Varandas &
C     J. Brandao Mol. Phys.,45,1982,857] without the 2-body correlation
C     energy.
         T = RR*(H2TB1 + RR*(H2TB2 + RR*(H2TB3 + RR*H2TB4)))
         AT = H2TA*EXP(-T)*RRI
         T = H2TB1 + RR*(2.0D0*H2TB2 + RR*(3.0D0*H2TB3 +
     *      RR*4.0D0*H2TB4))
         DAT = -AT*(T + RRI)
C     Add in triplet and subtract singlet 2-body correlation terms
         AT = AT + CORRT(I) - CORRS(I)
         DAT = DAT + DCORRT(I) - DCORRS(I)
C     Calculate effective diatomic triplet state curve.
         T = EXP(-AL3*RR)
         WE= (AL0 + RR*(AL1 + RR*AL2))*T
         DWE = (AL1 + 2.0D0*AL2*RR)*T - AL3*WE
C     Triplet energy
         ET = F*WE + OMF*AT
C     Contribution to coulomb term
         QSUM = QSUM + 0.5D0*(ES + ET)
C     Contribution to exchange term
         XJ(I) = 0.5D0*(ES - ET)
C     Derivatives of coulomb and exchange parts
         FTERM = WE - AT
         DO 20 J = 1,3
            T = 0.5D0*FTERM*DF(J)
            DQ(J) = DQ(J) + T
            DXJ(I,J) = -T
   20    CONTINUE
         T = F*DWE + OMF*DAT
         DQ(I) = DQ(I) + 0.5D0*(DESDR + T)
         DXJ(I,I) = DXJ(I,I) + 0.5D0*(DESDR - T)
   50 CONTINUE
C
      F1 = (XJ(1) - XJ(2))**2
      F2 = (XJ(2) - XJ(3))**2
      F3 = (XJ(3) - XJ(1))**2
      EXCH = SQRT(0.5D0*(F1 + F2 + F3))
      VLEPS = QSUM - EXCH
      VLEPS2 = QSUM + EXCH
      XTERM1 = 2.0D0*XJ(1) - XJ(2) - XJ(3)
      XTERM2 = 2.0D0*XJ(2) - XJ(1) - XJ(3)
      XTERM3 = 2.0D0*XJ(3) - XJ(1) - XJ(2)
      DO 80 I = 1,3
         TEXCH = EXCH
         IF(TEXCH.LE.0.0D0)TEXCH = 1.0D0
         XTERM = 0.5D0/TEXCH*(XTERM1*DXJ(1,I)  + XTERM2*DXJ(2,I)
     *      + XTERM3*DXJ(3,I))
         DLEP(I) = DQ(I) - XTERM
         DLEP2(I) = DQ(I) + XTERM
   80 CONTINUE
      RETURN
      END
      SUBROUTINE H3VA_HHH_DMBE (R1, R2, R3, VA, DVA)
C***********************************************************************
C     Calculates Va correction energy and its 1st derivative.
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION DVA(3)
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
      COMMON /H3COCM_HHH_DMBE/ PER, PER2, R12, R22, R32, QCOORD, DQ1, 
     $DQ2, DQ3,
     *   RHO, DRHO1, DRHO2, DRHO3, S, S2, DS1, DS2, DS3, CPHI3, DCPHI1,
     *   DCPHI2, DCPHI3
C
      T = ALPHA5*PER2
      EXPVA = EXP(-T*PER)
      IF (EXPVA .EQ. 0.0D0) THEN
         VA = 0.0D0
         DVA(1) = 0.0D0
         DVA(2) = 0.0D0
         DVA(3) = 0.0D0
      ELSE
         VA = 0.0D0
         DV = 0.0D0
         T3 = (R1-R2)*(R2-R3)*(R3-R1)
         T1 = T3*T3
         T2 = 1.0D0
         T4 = 2.0D0
         DO 1 J=1,4
            T2 = T2*T1
            VA = VA + XPAR(J)*T2
            DV = DV + T4*XPAR(J)*T3
            T3 = T3*T1
            T4 = T4 + 2.0D0
    1    CONTINUE
         VA = VA*EXPVA
         T1 = DV*EXPVA
         T2 = 3.0D0*T*VA
         DVA(1) = T1*(R2-R3)*(R3+R2-2.0D0*R1) - T2
         DVA(2) = T1*(R3-R1)*(R1+R3-2.0D0*R2) - T2
         DVA(3) = T1*(R1-R2)*(R1+R2-2.0D0*R3) - T2
      END IF
      RETURN
      END
      SUBROUTINE H3VII_HHH_DMBE (R1, R2, R3, E, DV2)
C***********************************************************************
C     Calculates VII correction energy and its 1st derivative.
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
      COMMON /H3COCM_HHH_DMBE/ PER, PER2, R12, R22, R32, QCOORD, DQ1, 
     $DQ2, DQ3,
     *   RHO, DRHO1, DRHO2, DRHO3, S, S2, DS1, DS2, DS3, CPHI3, DCPHI1,
     *   DCPHI2, DCPHI3
      DIMENSION DV2(3)
C
C   Compute B1 function.
      R1I = 1.0D0/R1
      R2I = 1.0D0/R2
      R3I = 1.0D0/R3
      COS1 = 0.5D0*R2I*R3I*(R12-R22-R32)
      COS2 = 0.5D0*R1I*R3I*(R22-R12-R32)
      COS3 = 0.5D0*R1I*R2I*(R32-R12-R22)
      WB = 1.0D0+COS1+COS2+COS3
      WB2 = WB*WB
C   WB derivatives
      WB1P = (R1*R3I - 1.0D0)*R2I - R3I - (COS2+COS3)*R1I
      WB2P = (R2*R3I - 1.0D0)*R1I - R3I - (COS1+COS3)*R2I
      WB3P = (R3*R2I - 1.0D0)*R1I - R2I - (COS1+COS2)*R3I
C   EB1 term
      EXP1 = EXP(-BETA1*PER)
      EXP3 = EXP(-BETA3*PER)
      EB1T = (XPAR(5) + XPAR(6)*PER)*EXP1
      EB3T = (XPAR(14) + XPAR(15)*PER2)*EXP3
      EB1 = WB*(EB1T + EB3T)
C   EB1 derivatives
      EB1PR = WB*(-BETA1*EB1T - BETA3*EB3T + XPAR(6)*EXP1
     *   + 2.0D0*PER*XPAR(15)*EXP3)
      EB1PWB = EB1T + EB3T
      EB1P1 = EB1PWB*WB1P + EB1PR
      EB1P2 = EB1PWB*WB2P + EB1PR
      EB1P3 = EB1PWB*WB3P + EB1PR
C   EB2 term
      T1 = BETA2*PER
      EXP2 = EXP(-T1*PER)
      EB2 = WB2*(XPAR(7) + WB*(XPAR(8) + WB*XPAR(9)))*EXP2
C   EB2 derivatives
      EB2PWB = WB*(2.0D0*XPAR(7) + WB*(3.0D0*XPAR(8) + WB*4.0D0*XPAR(9))
     V)*EXP2
      EB2PR = -2.0D0*T1*EB2
      EB2P1 = EB2PWB*WB1P + EB2PR
      EB2P2 = EB2PWB*WB2P + EB2PR
      EB2P3 = EB2PWB*WB3P + EB2PR
C   EB4 term
C      EB4A
      T2 = XPAR(10)*EXP1
      T3 = WB*XPAR(11)*EXP2
      EB4A = WB*(T2 + T3)
C      EB4A derivatives
      EB4APW = T2 + 2.0D0*T3
      EB4APR = -WB*(BETA1*T2 + 2.0D0*T1*T3)
      EB4AP1 = EB4APW*WB1P + EB4APR
      EB4AP2 = EB4APW*WB2P + EB4APR
      EB4AP3 = EB4APW*WB3P + EB4APR
C      EB4B
      T2 = XPAR(12)*EXP1
      T3 = XPAR(13)*EXP2
      EB4B = WB*(T2 + T3)
C      EB4B derivatives
      EB4BPW = T2 + T3
      EB4BPR = -WB*(BETA1*T2 + 2.0D0*T1*T3)
      EB4BP1 = EB4BPW*WB1P + EB4BPR
      EB4BP2 = EB4BPW*WB2P + EB4BPR
      EB4BP3 = EB4BPW*WB3P + EB4BPR
C
      DR12 = R1 - R2
      DR23 = R2 - R3
      DR31 = R3 - R1
      EQ = DR12*DR12 + DR23*DR23 + DR31*DR31
      EQP1 = 2.0D0*(DR12 - DR31)
      EQP2 = 2.0D0*(-DR12 + DR23)
      EQP3 = 2.0D0*(-DR23 + DR31)
      RI = R1I + R2I + R3I
      EB4 = EB4A*RI + EB4B*EQ
C   EB4 derivatives
      EB4P1 = EB4AP1*RI - EB4A/R12 + EB4BP1*EQ + EB4B*EQP1
      EB4P2 = EB4AP2*RI - EB4A/R22 + EB4BP2*EQ + EB4B*EQP2
      EB4P3 = EB4AP3*RI - EB4A/R32 + EB4BP3*EQ + EB4B*EQP3
C
      E = EB1 + EB2 + EB4
      DV2(1) = EB1P1 + EB2P1 + EB4P1
      DV2(2) = EB1P2 + EB2P2 + EB4P2
      DV2(3) = EB1P3 + EB2P3 + EB4P3
      RETURN
      END
      SUBROUTINE H3VIII_HHH_DMBE (R1, R2, R3, VIII, DV3)
C**********************************************************************
C     Calculates VIII correction energy and it 1st derivatives
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION DV3(3)
      COMMON /H3DMCM_HHH_DMBE/ ALPH2, ALPHA5, ALPH0, ALPH1, AL0, AL1, 
     $AL2, AL3,
     *   AZ2, BETA1, BETA2, BETA3, BET0, BET1, CD0, CD1, CHH(3), CK0(3),
     *   CK1(3), DAMPA(3), DAMPB(3), HFD, HFA1, HFA2, HFA3, HFGAM, H2RM,
     *   H2RMT, H2R0, H2TA, H2TB1, H2TB2, H2TB3, H2TB4, PI, RHOL, SQRT3,
     *   XPAR(15)
      COMMON /H3COCM_HHH_DMBE/ PER, PER2, R12, R22, R32, QCOORD, DQ1, 
     $DQ2, DQ3,
     *   RHO, DRHO1, DRHO2, DRHO3, S, S2, DS1, DS2, DS3, CPHI3, DCPHI1,
     *   DCPHI2, DCPHI3
C
      IF (S .EQ. 0.0D0) THEN
         VIII = 0.0D0
         DV3(1)=0.0D0
         DV3(2)=0.0D0
         DV3(3)=0.0D0
      ELSE
         T1 = RHO - RHOL
         T5 = ALPH2*T1
         EXPV3 = EXP(-T5*T1)
         IF (EXPV3 .EQ. 0.0D0) THEN
            VIII = 0.0D0
            DV3(1) = 0.0D0
            DV3(2) = 0.0D0
            DV3(3) = 0.0D0
         ELSE
            T1 = S2*CPHI3
            T2 = 1.0D0 + S*T1
            T3 = S2*T2
            T4 = CD0 + CD1*RHO
            VIII = T3*T4*EXPV3
            TDS   = (2.0D0*S*T2 + 3.0D0*S2*T1)*T4
            TDCPH = S2*S*S2*T4
            TDRHO = T3*(CD1 - 2.0D0*T5*T4)
            DV3(1) = (TDS*DS1 + TDCPH*DCPHI1 + TDRHO*DRHO1)*EXPV3
            DV3(2) = (TDS*DS2 + TDCPH*DCPHI2 + TDRHO*DRHO2)*EXPV3
            DV3(3) = (TDS*DS3 + TDCPH*DCPHI3 + TDRHO*DRHO3)*EXPV3
         END IF
      END IF
      RETURN
      END