C
C BELOW ARE RELEVANT CONSTANTS, VSM2 (THE MODIFIED SORBIE-MURRELL 
C SURFACE) AND VRADSM2 (DRIVER FOR OBTAINING POTENTIAL IN RADAU 
C COORDINATES).  I SEE THAT THE NOTATION IS A LITTLE INCONSISTENT, 
C BUT NOT BADLY.  THIS COMES FROM WORKING ON IT OVER TOO LONG A 
C PERIOD OF TIME.  UNITS ARE ANGSTROMS FOR LENGTH, AMU'S FOR MASS, 
C AND WAVE NUMBERS FOR ENERGY.  IF YOU FIND ANY MISTAKES, DON'T 
C LET ME KNOW. -- BRUCE JOHNSON
C**************************************************
C I HAVE CHANGED THE FIRST ARGUMENT TO SEND S3 (H-H DISTANCE)
C THROUGH DIRECTLY INSTEAD OF CTHETA=COS(THETA). - JDK 9/27/90
C**************************************************
C
CCCCCCFUNCTION VSM2(CTHETA,S1,S2)
      FUNCTION VSM2(S3,S1,S2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C RETURNS VALUE OF THE SORBIE-MURRELL WATER SURFACE FOR GIVEN VALUES
C OF THETA, THE INTERBOND ANGLE, AND S1 AND S2, THE TWO O-H BOND
C LENGTHS.
      DATA S3D0,DHHI,B1I,B2I,B3I,GI,BETA,DHHO,B1O,B2O,GO/
     1.7417D0,38311.8D0,3.87951D0,3.74086D0,3.25460D0,3.87951D0,
     215808.3D0,47238.2D0,1.47818D0,1.78015D0,3.9022D0/
      DATA C000I,C100I,C001I,C200I,C002I,C101I,C110I,C300I,C003I,C210I,
     1C102I,C201I,C111I,C400I,C004I,C310I,C220I,C301I,C202I,C103I,
     2C112I,C211I,C000O,C100O,C001O,C200O,C002O,C101O,C110O,C300O,C003O,
     3C210O,C102O,C201O,C111O,C400O,C004O,C310O,C220O,C301O,C202O,
     4C103O,C112O,C211O/1.D0,1.5402D0,4.6551D0,1.5748D0,-4.4369D0,
     517.995D0,-1.9861D0,6.0132D0,-7.5742D0,5.0266D0,11.391D0,
     622.279D0,23.089D0,-.5793D0,-11.478D0,.5220D0,-5.9434D0,
     759.011D0,-22.368D0,20.982D0,-18.712D0,93.800D0,1.,1.8676D0,
     8-0.3582D0,3.1773D0,-2.4096D0,3.1849D0,1.2402D0,6.8763D0,
     9-5.2016D0,8.1405D0,6.2925D0,-2.2399D0,-1.8101D0,-.3056D0,
     1-6.7552D0,.3061D0,-3.2029D0,8.8604D0,-15.311D0,8.1767D0,
     2-13.316D0,16.160D0/
      DATA S1E,S3E,S3P,AI,AO,GAM1,GAM3,GAMS/0.9572D0,1.5139D0,1.67067D0,
     1-7594.37D0,-13990.38D0,4.5348D0,2.D0,25.D0/
      DATA PI/3.14159265358979D0/
      DATA S1D0,S1P,AM12I,A0I,DOHO,A1O,A2O,A3O,F/0.9706D0,
     10.6758804D0,188.6254D0,-20756.86D0,37356.07D0,4.534735D0,
     25.03069D0,3.92551D0,4.534735D0/
C INITIALIZE REDMON-SCHATZ SWITCHING FACTOR
CCCCCCS3=DSQRT(S1*S1+S2*S2-2.D0*S1*S2*CTHETA)
      D3=S3-S3E
      DEL=S3-S3P
      IF(DEL.GT.1.6D0)THEN
      OMEGA=0.D0
      ELSE
      EE=DEXP(GAMS*(S3-S3P))
      OMEGA=PI/2.D0/(1.D0+EE*EE)
      ENDIF
      SSQ=DSIN(OMEGA)**2
      CSQ=1.D0-SSQ
      IF(DABS(GAM3*D3).GT.4.0D1)THEN
      FAC3=0.D0
      ELSE
      EE=DEXP(GAM3*D3)
      FAC3=2.D0/(1.D0+EE)
      ENDIF
C INITIALIZE FACTORS INVOLVING ONLY S1 AND S2 IF THESE ARE NEW
      D1=S1-S1E
      D2=S2-S1E
      D12=D1*D1
      D13=D12*D1
      D14=D13*D1
      D22=D2*D2
      D23=D22*D2
      D24=D23*D2
      IF(GAM1*D1/2.D0.GT.4.0D1)THEN
      FAC1=0.D0
      ELSE
      EE=DEXP(GAM1*D1/2.D0)
      FAC1=2.D0/(1.D0+EE*EE)
      ENDIF
      IF(GAM2*D2/2.D0.GT.4.0D1)THEN
      FAC2=0.D0
      ELSE
      EE=DEXP(GAM1*D2/2.D0)
      FAC2=2.D0/(1.D0+EE*EE)
      ENDIF
      F0I=C000I+C100I*(D1+D2)+C200I*(D12+D22)+C110I*D1*D2+C300I*(D13
     1+D23)+C210I*(D12*D2+D1*D22)+C400I*(D14+D24)+C310I*(D13*D2+D1*D23)
     2+C220I*D12*D22
      F1I=C001I+C101I*(D1+D2)+C201I*(D12+D22)+C111I*D1*D2+C301I*(D13
     1+D23)+C211I*(D12*D2+D1*D22)
      F2I=C002I+C102I*(D1+D2)+C202I*(D12+D22)+C112I*D1*D2
      F3I=C003I+C103I*(D1+D2)
      F4I=C004I
      F0O=C000O+C100O*(D1+D2)+C200O*(D12+D22)+C110O*D1*D2+C300O*(D13
     1+D23)+C210O*(D12*D2+D1*D22)+C400O*(D14+D24)+C310O*(D13*D2+D1*D23)
     2+C220O*D12*D22
      F1O=C001O+C101O*(D1+D2)+C201O*(D12+D22)+C111O*D1*D2+C301O*(D13
     1+D23)+C211O*(D12*D2+D1*D22)
      F2O=C002O+C102O*(D1+D2)+C202O*(D12+D22)+C112O*D1*D2
      F3O=C003O+C103O*(D1+D2)
      F4O=C004O
C COMBINE INNER AND OUTER POLYNOMIALS AND COEFFICIENTS FOR V3
      PREF=FAC1*FAC2*FAC3
      F0=AI*SSQ*F0I+AO*CSQ*F0O
      F1=AI*SSQ*F1I+AO*CSQ*F1O
      F2=AI*SSQ*F2I+AO*CSQ*F2O
      F3=AI*SSQ*F3I+AO*CSQ*F3O
      F4=AI*SSQ*F4I+AO*CSQ*F4O
      VSM2=F4*D3
      VSM2=(VSM2+F3)*D3
      VSM2=(VSM2+F2)*D3
      VSM2=(VSM2+F1)*D3
      VSM2=(VSM2+F0)*PREF
C CONSTRUCT VHH AND ADD TO VSM2
      T3=S3-S3D0
      VHHI=B3I*T3
      VHHI=(VHHI+B2I)*T3
      VHHI=(VHHI+B1I)*T3
      VHHI=BETA-DHHI*(1.D0+VHHI)*DEXP(-GI*T3)
      VHHO=B2O*T3
      VHHO=(VHHO+B1O)*T3
      VHHO=DHHO*(1.D0+VHHO)*DEXP(-GO*T3)
      VHH=SSQ*VHHI+CSQ*VHHO
      VSM2=VSM2+VHH
C CONSTRUCT THE O-H POTENTIALS FOR S1 AND S2 AND ADD TO VSM2
      S=S1
      N=0
10    CONTINUE
      IF(S.LT.S1P)THEN
      VOH=1.D0/S/S/S
      VOH=VOH*VOH*VOH*VOH
      VOH=VOH*AM12I+A0I
      ELSE
      T=S-S1D0
      VOH=T*A3O
      VOH=(VOH+A2O)*T
      VOH=(VOH+A1O)*T
      VOH=-DOHO*(1.D0+VOH)*DEXP(-F*T)
      ENDIF
      VSM2=VSM2+VOH
      N=N+1
      S=S2
      IF(N.EQ.1)GOTO 10
C ADD IN GROUND STATE ENERGY...
      VSM2=VSM2+81223.1015874772D0     
      RETURN
      END