SUBROUTINE H2H2
C
C     FARRAR-LEE H2-H2 MSV POTENTIAL PLUS
C       ANISOTROPIC TERMS USED BY ZARUR AND RABITZ.
C
C     ENTRY VINIT(I,RM,EPSIL)  INITIALIZES ROUTINE FOR LAM(I)
C                              SETS RM, EPSIL
C     ENTRY VSTAR (I,R,V)  GIVES POTENTIAL FOR LAM(I) AT R.
C     ENTRY VSTAR1(I,R,V) GIVES DERIVATIVE FOR LAM(I) AT R.
C     ENTRY VSTAR2(I,R,V) GIVES 2ND DERIV. FOR LAM(I) AT R.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE
      DIMENSION B(4),A(4,4)
      DIMENSION CONST(4)
C
C  DEFINE STATEMENT FUNCTIONS
      E(R)=EXP(BETA*(1D0-R))
      P4(C1,C2,C3,C4,R)=C1+R*(C2+R*(C3+R*C4))
      P3(C1,C2,C3,R)=C1+R*(C2+R*C3)
      P2(C1,C2,R)=C1+R*C2
C
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      ENTRY VINIT(I,RM,EPSIL)
C
C     SET CONSTANT FACTORS . . .
C
      RM=3.49D0
      EPSIL=24.17D0
      RMSAVE=RM
      CONST(1)=(4D0*3.1415926D0)**( 1.5D0)
      CONST(2)=(.14D0/5D0)*CONST(1)
      CONST(3)=CONST(2)
      CONST(4)=54836D0 / EPSIL
C
      IF (I.EQ.4)  GO TO 1400
      IF (I.LT.1  .OR. I.GT.4)  GO TO 9999
C
      WRITE(6,601)
  601 FORMAT('0',10X,'FARRAR-LEE MSV H2-H2 POTENTIAL /     RM=3.49,',
     1        ' EPSIL=24.17   SET INTERNALLY.')
      IF (I.GT.1)  WRITE(6,606)
  606 FORMAT(16X,'ANISOTROPIC TERM IS .14/5.0 TIMES ISOTROPIC TERM')
C
C     THE FOLLOWING PARAMETERS MUST BE SPECIFIED INTERNALLY.
C        R1,R2 IN RM .   C6,C8 IN (1/CM )/ANG**N        BETA
C
C     SCALE PARAMETERS WITH RM, EPS.
C
      R1=1.0754D0
      R2=1.4D0
      C6=( 57913.D0 /EPSIL)/  RM**6
      C8=( 156263.D0/EPSIL)/  RM**8
      BETA=6.5D0
C
      WRITE(6,602)  R1,BETA,R2,C6,C8
  602 FORMAT(20X,'FOR R LESS THAN',F8.4,', MORSE BETA =',F9.4/
     1  20X,'FOR R GREATER THAN',F8.4,', VAN DER WAAL C6, C8 =',
     2      2F10.5,' (EPSIL/RM**N).')
C
C     SET-UP SPLINE COEFFS. B(J,I), J=1,4  FOR LAM(I).
C     ** MAR 1979 CHANGED TO INCORPORATE PRECOMPUTED SPLINE COEFFS.
C
      B(1)=-.6829 4444 3121 D1
      B(2)= .7211 3009 9412 D1
      B(3)=-.7689 9426 2915 D0
      B(4)=-.7125 2209 5040 D0
C
 1100 WRITE(6,605)  (B(J  ),J=1,4)
  605 FORMAT(20X,'CUBIC SPLINE COEFFICIENTS, TO INCREASING POWERS OF R,
     1ARE AS FOLLOWS'/25X,4D16.8)
      RETURN
C
 1400 WRITE(6,607)
  607     FORMAT('0',10X,'QUADRUPOLE-QUADRUPOLE LONG-RANGE W/ Q=.662 BUC
     1KINGHAM.')
      AC=.0687D0
      RETURN
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      ENTRY VSTAR (I,R,V)
      IF (I.EQ.4)  GO TO 2400
      IF (R.LE.R1)  GO TO 2100
      IF (R.GE.R2)  GO TO 2200
      V=P4(B(1),B(2),B(3),B(4),R)
      GO TO 5000
 2100 T1=E(R)
      V=T1*(T1-2D0)
      GO TO 5000
 2200 RSQ=1D0/(R*R)
      R6=RSQ*RSQ*RSQ
      V=-       R6*(C6+C8*RSQ)
      GO TO 5000                                                          99.
 2400 RR=R*RMSAVE
      R12=R**(-12)
      RRT=RR+AC*R12
      V=RRT**(-5)
      GO TO 5000
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      ENTRY VSTAR1(I,R,V)
      IF (I.EQ.4)  GO TO 3400
      IF (R.LE.R1)  GO TO 3100
      IF (R.GE.R2)  GO TO 3200
      V=P3(B(2),2D0*B(3),3D0*B(4),R)
      GO TO 5000
 3100 T1=E(R)
      V=-2D0*BETA*       T1*(T1-1D0)
      GO TO 5000
 3200 RSQ=1D0/(R*R)
      R6=RSQ*RSQ*RSQ/R
      V=       R6*(6D0*C6+8D0*C8*RSQ)
      GO TO 5000
 3400 RR=R*RMSAVE
      R12=R**(-12)
      RRT=RR+AC*R12
      V=-5D0*RRT**(-6)*(1D0-12D0*AC*R12/RR)
      GO TO 5000
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      ENTRY VSTAR2(I,R,V)
      IF (I.EQ.4)  GO TO 4400
      IF (R.LE.R1)  GO TO 4100
      IF (R.GE.R2)  GO TO 4200
      V=P2(2D0*B(3),6D0*B(4),R)
      GO TO 5000
 4100 T1=E(R)
      V=T1*BETA*BETA*       (4D0*T1-2D0)
      GO TO 5000
 4200 RSQ=1D0/(R*R)
      R6=RSQ*RSQ
      R6=R6*R6
      V=-       R6*(42D0*C6+56D0*C8*RSQ)
      GO TO 5000
 4400 RR=R*RMSAVE
      R12=R**(-12)
      RRT=RR+AC*R12
      XXX=1D0-12D0*AC*R12/RR
      V=-5D0*RRT**(-6)*(156D0*AC*R12/(RR*RR)-6D0*XXX*XXX/RRT)
      GO TO 5000
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     UNIFIED RETURN POINT
 5000 V=V*CONST(I)
      RETURN
C
C
 9999 WRITE(6,699)  I
  699 FORMAT('0 * * * ERROR.  POTENTIAL NOT DEFINED FOR SYMMETRY=',I10)
      STOP
      END
      SUBROUTINE VRTP(IDERIV,R,VV)
C
C     FARRAR-LEE H2-H2 MSV POTENTIAL PLUS
C       ANISOTROPIC TERMS USED BY ZARUR AND RABITZ.
C     COMMON /ANGLES/ FOR VERSION 14
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE
      DIMENSION B(4),A(4,4)
      DIMENSION CONST(4),VI(4),FACT(4)
C
      COMMON/ANGLES/COSANG(7),FACTOR,IHOMO,ICNSYM,IHOMO2,ICNSY2
C
C  DEFINE STATEMENT FUNCTIONS
      E(R)=EXP(BETA*(1D0-R))
      P4(C1,C2,C3,C4,R)=C1+R*(C2+R*(C3+R*C4))
      P3(C1,C2,C3,R)=C1+R*(C2+R*C3)
      P2(C1,C2,R)=C1+R*C2
C
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      IF (IDERIV.EQ.-1) THEN
        RM=3.49D0
        RMSAVE=RM
        EPSIL=24.17D0
C       SET UP CONSTANTS
        CONST(1)=(4D0*3.1415926D0)**( 1.5D0)
        CONST(2)=(.14D0/5D0)*CONST(1)
        CONST(3)=CONST(2)
        CONST(4)=54836D0 / EPSIL
        R1=1.0754D0
        R2=1.4D0
        C6=( 57913.D0 /EPSIL)/  RM**6
        C8=( 156263.D0/EPSIL)/  RM**8
        BETA=6.5D0
        AC=.0687D0
C       SET-UP SPLINE COEFFS. B(J,I), J=1,4  FOR LAM(I).
        B(1)=-.6829 4444 3121 D1
        B(2)= .7211 3009 9412 D1
        B(3)=-.7689 9426 2915 D0
        B(4)=-.7125 2209 5040 D0
        WRITE(6,601) RM,EPSIL
  601   FORMAT(10X,'FARRAR-LEE MSV H2-H2 POTENTIAL'/10X,'RM=',F8.3/
     1         10X,'EPSIL=',F12.2)
        WRITE(6,606)
  606   FORMAT(10X,'ANISOTROPIC TERM IS .14/5.0 TIMES ISOTROPIC TERM')
        WRITE(6,602)  R1,BETA,R2,C6,C8
  602   FORMAT(10X,'FOR R LESS THAN',F8.4,', MORSE BETA =',F9.4/
     1    10X,'FOR R GREATER THAN',F8.4,', VAN DER WAAL C6, C8 =',
     2      2F10.5,' (EPSIL/RM**N).')
        WRITE(6,607)
  607   FORMAT(10X,'QUADRUPOLE-QUADRUPOLE LONG-RANGE W/ Q=.662 BU',
     1         'CKINGHAM.')
        WRITE(6,605)  (B(J),J=1,4)
  605   FORMAT(10X,'CUBIC SPLINE COEFFICIENTS, TO INCREASING ',
     1           'POWERS OF R, ARE AS FOLLOWS'/15X,4D16.8)
        PI=ACOS(-1.D0)
        PIF=4.D0*PI
        FACT(1)=1.D0/SQRT(PIF*PIF*PIF)
        FACT(2)=FACT(1)*2.5D0
        FACT(3)=FACT(1)
        FACT(4)=FACT(1)*(45.D0/4.D0)*SQRT(1.D0/7.D1)
C
        IHOMO=2
        ICNSYM=2
C       RETURN RM EPSIL
        R=RM
        VV=EPSIL
        RETURN
C
      ELSEIF (IDERIV.EQ.0) THEN
C
        IF (R.LE.R1)  GO TO 2100
        IF (R.GE.R2)  GO TO 2200
        V=P4(B(1),B(2),B(3),B(4),R)
        GO TO 5000
 2100   T1=E(R)
        V=T1*(T1-2D0)
        GO TO 5000
 2200   RSQ=1D0/(R*R)
        R6=RSQ*RSQ*RSQ
        V=-R6*(C6+C8*RSQ)
 5000   DO 5001 I=1,3
 5001   VI(I)=V*CONST(I)
C       QUADRUPOLE - QUADRUPOLE TERM
        RR=R*RMSAVE
        R12=R**(-12)
        RRT=RR+AC*R12
        V=RRT**(-5)
        VI(4)=V*CONST(4)
C       CT1=COSANG(1)
C       ST1=SQRT(1.D0-CT1*CT1)
C       CT2=COSANG(2)
C       ST2=SQRT(1.D0-CT2*CT2)
C       P21=3.D0*CT1*CT1-1.D0
C       P22=3.D0*CT2*CT2-1.D0
C       ADD THE PIECES WITH COSANG() VALUES - YRR GIVES ROTOR FUNCTIONS
        VV=VI(1)*YRR(0,0,0,COSANG(1),COSANG(2),COSANG(3))
     1    +VI(2)*YRR(2,0,2,COSANG(1),COSANG(2),COSANG(3))
     2    +VI(3)*YRR(0,2,2,COSANG(1),COSANG(2),COSANG(3))
     1    +VI(4)*YRR(2,2,4,COSANG(1),COSANG(2),COSANG(3))
        RETURN
 
      ELSE
C
        WRITE(6,*) ' *** VRTP (H2-H2). NO SUPPORT FOR IDERIV',IDERIV
        STOP
      ENDIF
      END