CCL Home Page
Up Directory CCL spline
      SUBROUTINE POTENL(ICNTRL, MXLMB, MPOTL, LAM, R, P, ITYP)
C
C     THIS IS A REPLACEMENT ROUTINE FOR SUBROUTINE POTENL IN MOLSCAT.
C     IT IS DESIGNED TO HANDLE CASES WHERE THE POTENTIAL ANGULAR 
C       EXPANSION TERMS, V_SUB_ISYM, ARE GIVEN BY A TABLE OF VALUES
C       AT RADIAL DISTANCES R(I)=RMIN+(I-1)*DR, I=1,NPT.      
C     N.B. THE SAME, EQUALLY SPACED RADIAL GRID MUST BE AVAILABLE
C       FOR ALL THE ANGULAR SYMETRIES.
C     DATA ARE READ FROM UNIT IIN (DEFAULT=12) USING FORMAT 
C       STATEMENTS 501, 502, 503.  SEE READ AND FORMAT STATEMENTS
C       FOR A DESCRIPTION OF THE INPUT.
C     CONTINUOUS RADIAL FUNCTIONS ARE OBTAINED FROM THE INPUT POINTS
C       BY CUBIC SPLINE INTERPOLATION, EXPONENTIAL EXTRAPOLATION AT
C       SMALL DISTANCES, AND INVERSE POWER EXTRAP AT LARGE DISTANCES.
C       (S. GREEN, J CHEM PHYS 67, 715 (1977).
C
C     **************************************************************
C     ** THIS DECK COMPATIBLE WITH MOLSCAT VERSION 14 (FEB 94)    **
C     ** USES IV INDEXING W/ ITYPE=2                              **
C     ** USES UPPER /MEMORY/ FOR STORAGE. IT IS POSSIBLE TO       **
C     **   DESTROY LAM() BY OVERWRITING W/ Y()..SCR() BUT THIS    **
C     **   SHOULD BE CAUGHT BY CHKSTR AFTER RETURN FROM POTENL    **
C     ** NB MANY STANDARD &POTL NAMELIST ITEMS HAVE BEEN REMOVED  **
C     **   DESCRIPTION OF VARIABLES FOLLOWS NAMELIST STATEMENT    **
C     **************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE IXMEM,RMINSV,DXSV,NPSV,RLSV,RUSV
C
      PARAMETER (IIN=12)
      PARAMETER (MXLIN=5)
C
      DIMENSION P(MXLMB), LAM(MXLMB)
      DIMENSION LAMIN(MXLIN)
      DIMENSION SP(3)
      INTEGER CFLAG
      CHARACTER*8 QNAME(10), QTYPE(10)
      CHARACTER*8 LABEL(10)
      EQUIVALENCE (MXLAM,MXSYM)
      COMMON/NPOT/NPTL
C
      COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1)
C
      DATA MXSKPR/50/
      DATA QTYPE/'LAMBDA =','ABS(MU)=','   MU = ','   L1 = ',
     1           '   L2 = ','    L = ','    V = ','V-PRIME=',
     2           '   J =  ','J-PRIME='/
C
      NAMELIST/POTL/RM,EPSIL, MXLAM,MXSYM,       NPOTL,CFLAG,
     1   LAMMAX,MXV,MXDELV, IPRINT
C
C     RM     - LENGTH SCALING FACTOR
C     EPSIL  - ENERGY SCALING FACTOR
C     MXLAM  - NUMBER OF POTENTIAL TERMS RETURNED
C     MXSYM  - SYNONYM FOR MXLAM (EQUIVALENCED VARIALBES)
C     NPOTL  - NUMBER OF ELEMENTS OF VL ARRAY (IN BASE) WHICH
C              CORRESPOND TO EACH ELEMENT OF P ARRAY
C     CFLAG  - FLAG FOR SCALING POTENTIAL FOR ITYPE = 5 OR 6
C     IPRINT - CONTROLS OUTPUT OF POTENTIAL INFORMATION
C              =0 (MINIMAL), =1 (SUMMARY), =2 (FULL)
C     ------ VARIABLES TO CONTROL RETAINING RECORDS FROM FILE(IIN)
C     LAMMAX  - HIGHEST LEGENDRE TERM TO KEEP
C     MXV    - HIGHEST VIBRATIONAL LEVEL TO KEEP (ITYPE=2)
C     MXDELV - HIGHEST DELTA-V TO KEEP (ITYPE=2)
C
      IF (ICNTRL.GE.0.AND.ICNTRL.LE.2) GOTO 1000
      IF (ICNTRL.EQ.-1) GOTO 2000
      WRITE(6,631) ICNTRL,R
  631 FORMAT('0 *** ERROR IN POTENL, ICNTRL =',I6,'  R =',E16.8)
      STOP
C
C     EVALUATE V(R)
C
 1000 NSTOR=2*NPSV+4
      IXSC=IXMEM-NPSV
      IF (R.LT.RLSV) GO TO 5911
C     IF (R.GT.RUSV) GO TO 5912 - BUG FIXED 3 MAR 95, SG
      IF (R.GE.RUSV) GO TO 5912
C
C     BELOW FOR CASE THAT POTENTIAL IS INTERPOLATED.
      XMIN=RMINSV
      DR=DXSV
      N=NPSV
C     FIND I SUCH THAT R IS IN INTERVAL X(I) TO X(I+1).
      Q=(R-XMIN)/DR
      IF (Q.LT.0D0)  GO TO 9999
      I=Q
      I=I+1
      IF (I.GE.N)  GO TO 9999
      HT1=R-XMIN-(I-1)*DR
      HT2=HT1-DR
      Z=HT1*HT2
C     LOOP OVER SYMMETRIES
      DO 1111 J=1,MXLMB
      IXY=IXSC-(NPSV+4)
      SSS=(X(IXSC+I+1)-X(IXSC+I))/DR
C     CALCULATE 2ND DERIVATIVE AT T.
      SP(3)=X(IXSC+I)+HT1*SSS
C     CALCULATE FUNCTION AT T.
      DELSQS=(X(IXSC+I)+X(IXSC+I+1)+SP(3))/6D0
      DELY=(X(IXY+I+1)-X(IXY+I))/DR
      SP(1)=X(IXY+I)+HT1*DELY   +Z*DELSQS
C     FIND 1ST DERIVATIVE AT T.
      SP(2)=DELY+(HT1+HT2)*DELSQS+Z*SSS/6D0
      P(J)=SP(ICNTRL+1)
 1111 IXSC=IXSC-NSTOR
      RETURN
 9999 WRITE(6,692)    R
  692 FORMAT('0 *** POTENL(SPLINE) ERROR INTERP REQUESTED FOR '
     1      ,'ILLEGAL R =',D16.8)
      STOP
C
C     EXTRAPOLATION AT LEFT
 5911 DO 5997 I=1,MXLMB
      TERM=0.D0
      AA=X(IXSC-3)
      B=X(IXSC-2)
      IF (AA   .EQ.0.D0 .OR. B   .EQ.0.D0) GO TO 5897
      TERM=AA   *EXP(-B   *R)
      IF (ICNTRL.GT.0) TERM=TERM*(-B)**ICNTRL
 5897 P(I)=TERM
 5997 IXSC=IXSC-NSTOR
      RETURN
C
C     EXTRAPOLATE AT RIGHT
 5912 DO 5996 I=1,MXLMB
      C=X(IXSC-1)
      D=X(IXSC)
      TERM=0.D0
      IF (C   .EQ.0.D0 .OR. D   .EQ.0.D0) GO TO 5896
      TERM=C   *R**(-D   )
      IF(ICNTRL.EQ.1) TERM=-D   *TERM/R
      IF(ICNTRL.EQ.2) TERM=D   *(D   +1.D0)*TERM/(R*R)
 5896 P(I)=TERM
 5996 IXSC=IXSC-NSTOR
      RETURN
C
C     POTENTIAL INITIALISATION
C
 2000 WRITE(6,630) IIN
  630 FORMAT('0 POTENL(SPLINE) ROUTINE (FEB 94).'/
     1 '0 TABULATED V(ISYM) WILL BE INPUT FROM UNIT',I4)
C
      RM=0.D0
      EPSIL=0.D0
      NPOTL=-1
      MXLAM=-1
      LAMMAX=-1
      MXV=-1
      MXDELV=-1
      CFLAG=0
      IPRINT=1
      READ(5,POTL)
C
C     REPORT &POTL INPUT AND CHECK CONSISTENCY W/ POTENL(SPLINE)
C
      ITYPE=ITYP-10*(ITYP/10)
      WRITE(6,634) IPRINT
  634 FORMAT('0 PRINT LEVEL FROM &POTL IPRINT =',I3)
      IF(NPOTL.NE.-1) THEN
        WRITE(6,*) ' *** POTENL(SPLINE) IGNORING INPUT NPOTL',NPOTL
      ENDIF
      IF (RM.NE.0.D0.OR.EPSIL.NE.0.D0) THEN
        WRITE(6,*) ' *** POTENL(SPLINE) INPUT RM, EPSIL =',RM,EPSIL
        WRITE(6,*) '     WILL BE IGNORED IN FAVOR OF VALUES ON (IIN).'
        WRITE(6,*) ' *** BE SURE THIS IS CORRECT.'
      ENDIF
      IF (MXLAM.GT.0) THEN
        WRITE(6,*) ' *** POTENL(SPLINE) INPUT MXLAM =',MXLAM
        WRITE(6,*) '     MAY LIMIT NUMBER OF SYMMETRIES READ'
      ENDIF
      IF (LAMMAX.GE.0) THEN
        WRITE(6,*) ' *** POTENL(SPLINE) INPUT LAMMAX =',LAMMAX
        WRITE(6,*) '     MAY LIMIT NUMBER OF SYMMETRIES RETAINED'
      ENDIF
      IF (MXV.GE.0.AND.(ITYPE.EQ.2.OR.ITYPE.EQ.7)) THEN
        WRITE(6,*) ' *** POTENL(SPLINE) INPUT MXV =',MXV
        WRITE(6,*) '     MAY LIMIT NUMBER OF SYMMETRIES RETAINED'
      ENDIF
      IF (MXDELV.GE.0.AND.(ITYPE.EQ.2.OR.ITYPE.EQ.7)) THEN
        WRITE(6,*) ' *** POTENL(SPLINE) INPUT MXDELV =',MXDELV
        WRITE(6,*) '     MAY LIMIT NUMBER OF SYMMETRIES RETAINED'
      ENDIF
C
C     GET AND CHECK HEADER INFORMATION FROM UNIT(IIN)
C       MXLTP IS NO. SYMMETRIES ON TAPE; NQPLTP IS NO. LABELS/SYMMETRY
C       EACH SYMMETRY DESCRIBED BY NP VALUES FROM XMIN IN STEPS OF DR
C       RM, EPSIL ARE LENGTH AND ENERGY SCALING VALUES FOR MOLSCAT
C
      READ(IIN,501,END=9100) LABEL,MXLTP,NQPLTP,NP,XMIN,DR,RM,EPSIL
  501 FORMAT(10A8/3I5,4F10.4)
      WRITE(6,632) LABEL
  632 FORMAT(/'  LABEL FROM POTENTIAL FILE --'/1X,10A8)
      NPSV=NP
      RMINSV=XMIN
      DXSV=DR
      RLSV=XMIN
      RUSV=XMIN+(NP-1)*DR
      IF (MXLAM.GE.1.AND.MXLAM.LT.MXLTP) THEN
        WRITE(6,*) ' *** POTENL(SPLINE).  INPUT MXLAM =',MXLAM
        WRITE(6,*) '     IS LESS THAN VALUE FROM TAPE',MXLTP
        WRITE(6,*) '     AND WILL LIMIT NUMBER OF SYMMETRIES INPUT'
        MXLTP=MXLAM
      ENDIF
      IF (IPRINT.LT.2) WRITE(6,633) NPSV,RMINSV,DXSV,RUSV
  633 FORMAT(/'  INTERPOLATION ON',I4,' POINTS FROM',F8.3,
     1       ' IN STEPS OF',F7.3,' TO',F8.3)
      IF (NPSV.LT.2) THEN
        WRITE(6,*) ' *** POTENL(SPLINE).  ERROR.  TOO FEW POINTS'
     1   ,' FOR SPLINE',NPSV
        STOP
      ENDIF
C
C     ATTEMPT TO PROCESS ITYPE AND POTENTIAL DESCRIPTION NUMBERS
C
      IF (IPRINT.GE.1) WRITE(6,639)
  639 FORMAT(/'  ANGULAR DEPENDENCE OF POTENTIAL EXPANDED IN TERMS OF')
      IF(ITYPE.EQ.1) GOTO 2001
      IF(ITYPE.EQ.2) GOTO 2002
      IF(ITYPE.EQ.3) GOTO 2003
      IF(ITYPE.EQ.5) GOTO 2005
      IF(ITYPE.EQ.6) GOTO 2005
      IF(ITYPE.EQ.7) GOTO 2002
C
      WRITE(6,640) ITYPE
  640 FORMAT('  *** POTENL(SPLINE). ITYPE =',I4,' IS NOT SUPPORTED.')
      STOP
C
 2001 NQPL=1
      QNAME(1)=QTYPE(1)
      IF (IPRINT.GE.1) WRITE(6,641)
  641 FORMAT('   LEGENDRE POLYNOMIALS, P(LAMBDA).')
      GOTO 2100
C
 2002 NQPL=3
      QNAME(1)=QTYPE(1)
      QNAME(2)=QTYPE(7)
      QNAME(3)=QTYPE(8)
      IV1=2
      IV2=3
      IF (IPRINT.GE.1) THEN
       WRITE(6,641)
       WRITE(6,642)
      ENDIF
  642 FORMAT('  INTEGRATED OVER DIATOM VIBRATIONAL FUNCTIONS')
      IF(ITYPE.NE.7) GOTO 2100
      NQPL=5
      QNAME(3)=QTYPE(9)
      QNAME(4)=QTYPE(8)
      QNAME(5)=QTYPE(10)
      IV2=4
      IF (IPRINT.GE.1) WRITE(6,643)
  643 FORMAT('  FOR EACH PAIR OF V,J LEVELS')
      GOTO 2100
C
 2003 NQPL=3
      QNAME(1)=QTYPE(4)
      QNAME(2)=QTYPE(5)
      QNAME(3)=QTYPE(6)
      IF (IPRINT.GE.1) WRITE(6,644)
  644 FORMAT('  CONTRACTED NORMALISED SPHERICAL HARMONICS, SUM',
     1   '(M1,M2,M) C(L1,M1,L2,M2,L,M) Y(L1,M1) Y(L2,M2) Y(L,M)'/
     2   '  SEE RABITZ, J. CHEM. PHYS. 57, 1718 (1972)')
      GOTO 2100
C
 2005 NQPL=2
      QNAME(1)=QTYPE(1)
      QNAME(2)=QTYPE(2)
      IF (IPRINT.GE.1) WRITE(6,645)
  645 FORMAT('  NORMALISED SPHERICAL HARMONICS: (Y(LAM,MU) + ',
     1       '(-)**MU Y(LAM,-MU)) / (1+DELTA(MU,0))')
      IF(CFLAG.EQ.1) WRITE(6,646)
  646 FORMAT('  *** POTENL(SPLINE) CFLAT=1 REQUEST TO SCALE POTENTIAL',
     1   ' IGNORED.'/'  *** MAKE SURE THIS IS OKEY.')
      GOTO 2100
C
 2100 IF (NQPLTP.NE.NQPL) THEN
        WRITE(6,*) ' *** POTENL(SPLINE) - POSSIBLE ERROR'
        WRITE(6,*) '     NQPL APPROPRIATE TO ITYPE',ITYPE
        WRITE(6,*) '     NOT EQUAL TO VALUE ON TAPE',NQPLTP
        WRITE(6,*) '     LATTER WILL BE USED'
        NQPL=NQPLTP
        IF (NQPL.GT.MXLIN) THEN
          WRITE(6,*) ' *** ERROR. EXCEEDS DIMENSION MXLIN',MXLIN
          STOP
        ENDIF
      ENDIF
C
C     GET INFORMATION FOR EACH SYMMETRY FROM UNIT(IIN)
      IQ=0
      MXLAM=0
      NSKIP=0
C     PREPARE TO USE /MEMORY/...,X  -- ALLOCATE SPACE FOR X(NPSV)
      IXMEM=MX
      MX=MX-NPSV
      DO 9000 I=1,MXLTP
      READ(IIN,502,END=9100) (LAMIN(IX),IX=1,NQPL)
  502 FORMAT(10I5)
C     SEE IF THIS SET MEETS ALL REQUIRMENTS OF LAMMAX,MXV,MXDELV
      IF (LAMMAX.GE.0.AND.LAMIN(1).GT.LAMMAX) THEN
      NSKIP=NSKIP+1
        IF (IPRINT.GE.1.AND.NSKIP.LE.MXSKPR) THEN
          WRITE(6,*) ' *** POTENL(SPLINE). SKIPPING INPUT SET BECAUSE',
     1               ' OF LAMMAX',LAMMAX
          WRITE(6,*) (LAMIN(IX),IX=1,NQPL)
        ENDIF
        GO TO 9001
      ENDIF
      IF (ITYPE.EQ.2.OR.ITYPE.EQ.7) THEN
        IF (MXV.GE.0.AND.(LAMIN(IV1).GT.MXV.OR.LAMIN(IV2).GT.MXV)) THEN
        NSKIP=NSKIP+1
          IF (IPRINT.GE.1.AND.NSKIP.LE.MXSKPR) THEN
            WRITE(6,*) ' *** POTENL(SPLINE). SKIPPING INPUT SET',
     1                 ' BECAUSE OF MXV',MXV
            WRITE(6,*) (LAMIN(IX),IX=1,NQPL)
          ENDIF
        GO TO 9001
        ENDIF
        IF (MXDELV.GE.0.AND.ABS(LAMIN(IV1)-LAMIN(IV2)).GT.MXDELV) THEN
        NSKIP=NSKIP+1
          IF (IPRINT.GE.1.AND.NSKIP.LE.MXSKPR) THEN
            WRITE(6,*) ' *** POTENL(SPLINE). SKIPPING INPUT SET',
     1                 ' BECAUSE OF MXDELV',MXDELV
            WRITE(6,*) (LAMIN(IX),IX=1,NQPL)
          ENDIF
        GO TO 9001
        ENDIF
      ENDIF
C
C     IF WE REACH HERE, RETAIN THE INPUT SYMMETRY.  CHECK STORAGE
      MXLAM=MXLAM+1
      IF (MXLAM*NQPL.GT.MXLMB) THEN
        WRITE(6,*) ' *** POTENL(SPLINE). INADEQUATE STORAGE IN EXTERNAL'
     1            ,' LAM ARRAY'
        WRITE(6,*) '     MXLAM, NQPL, MXLMB =',MXLAM,NQPL,MXLMB
        STOP
      ENDIF
C     STORAGE FOR Y(NPSV),AA,B,C,D,SCR(NPSV)
      NSTOR=2*NPSV+4
      MX=MX-NSTOR
      IF (MX.LT.IXNEXT) THEN
        WRITE(6,*) ' *** POTENL(SPLINE) INSUFFICIENT SPACE IN /MEMORY/'
        WRITE(6,*) '     TRYING TO PROCESS SYMMETRY',MXLAM
        WRITE(6,*) '     ORIGINAL/CURRENT MX, IXNEXT =',IXMEM,MX,IXNEXT
        STOP
      ENDIF
      IXX=MX
      IXY=MX+NPSV
      IXSC=IXY+NPSV+4
C
C     PUT LAMIN() VALUES INTO LAM() AND SET UP, IE, INPUT Y(,MXLAM)
      DO 9010 J=1,NQPL
 9010 LAM(IQ+J)=LAMIN(J)
      IF (IPRINT.GE.1) THEN
        WRITE(6,651) MXLAM
  651   FORMAT('0 INTERACTION POTENTIAL FOR SYMMETRY TYPE NUMBER',I3)
        WRITE(6,652) (QNAME(J),LAM(IQ+J),J=1,NQPL)
  652   FORMAT('  WHICH HAS ',6(A8,I3,3X))
      ENDIF
      READ(IIN,503,END=9100) (X(IXY+J),J=1,NPSV)
  503 FORMAT(4D20.12)
      DO 1011 J=1,NPSV
 1011 X(IXX+J)=RMINSV+(J-1)*DXSV
      IF (IPRINT.GE.2) THEN
        WRITE(6,654)
  654   FORMAT(1X)
        WRITE(6,601) NPSV,(X(IXX+J),X(IXY+J),J=1,NPSV)
  601   FORMAT('0',5X,'POTENTIAL WILL BE INTERPOLATED FROM FOLLOWING',
     1   I5,'  POINTS'/ ( 9X,4(F7.4,F15.6,8X)))
      ENDIF
C
C     GET EXTRAPOLATION TO LEFT AS Y =A(I)*DEXP(-B(I)*R)
      Y1=X(IXY+1)
      Y2=X(IXY+2)
      X1=X(IXX+1)
      X2=X(IXX+2)
      IF (Y1*Y2.GT.0.)  GO TO 1300
      IF (IPRINT.GE.1) WRITE(6,697)
  697 FORMAT('0 * * * WARNING.  SHORT RANGE CANNOT BE FIT AS EXPONENTIAL
     &.  SET TO ZERO.')
      AA       =0.
      B       =0.
      SCR1=0.D0
      GO TO 1399
 1300 B       =DLOG(Y1/Y2)/(X2-X1)
      AA       =Y1*DEXP(B*X1)
      SCR1=AA*B*B*DEXP(-B*X1)
 1399 IF (IPRINT.GE.2) WRITE(6,602)  RLSV,AA       ,B
  602 FORMAT('0  FOR R LESS THAN',F7.3,', V =',1P,D12.4,' * DEXP( -',
     1   0P,F8.3,' * R ).')
C
C     AND THEN EXTRAPOLATION TO RIGHT AS C*R**(-D)
      YN=X(IXY+NPSV)
      YNM1=X(IXY+NPSV-1)
      XN=X(IXX+NPSV)
      XNM1=X(IXX+NPSV-1)
      IF (YNM1*YN.LE.0D0) GOTO 1299
      IF (YNM1/YN.GT.1D0)  GO TO 1200
 1299 IF (IPRINT.GE.1) WRITE(6,698)
  698 FORMAT('0 * * * WARNING.  LONG-RANGE CANNOT BE FIT AS DECREASING E
     1XPONENTIAL.  SET TO ZERO.')
      D       =0D0
      C       =0D0
      SCRN=0.D0
      X(IXY+NPSV)=0.D0
      GO TO 1100
 1200 D       =
     1  DLOG(YNM1/YN)/DLOG(XN/XNM1)
C     ELIMINATE LONG=RANGE THAT DECREASES EXTREMELY FAST
      IF (D       .GT.8.D1)  GO TO 1299
      C       =YN*XN**D
      SCRN           =D       *(D       +1.D0)*
     1                C       *XN**(-D       -2.D0)
C     WEED OUT POSSIBLY SPURIOUS FITS WHICH FALL OFF TOO SLOWLY.
      IF (D       .GT.4.D0)  GO TO 1100
      IF (IPRINT.GE.1) WRITE(6,686)  D
  686 FORMAT('0 * * * WARNING.  INVERSE POWER =',D14.6,'  IS UNPHYSICAL.
     1  LONG RANGE FORCED TO ZERO UNLESS SPLNEQ CHANGED.')
      GO TO 1299
 1100 IF (IPRINT.GE.2) WRITE(6,603) RUSV,C       ,D
  603 FORMAT('0  FOR R GREATER THAN',F7.3,', V =',1P,D12.4,' / R**(',
     1  0P,F8.3,' ).')
C     BUT AA,B,C,D,SCR1,SCRN INTO APPROPRIATE STORAGE
      X(IXSC-3)=AA
      X(IXSC-2)=B
      X(IXSC-1)=C
      X(IXSC)=D
      X(IXSC+1)=SCR1
      X(IXSC+NPSV)=SCRN
C     INITIALIZE SCR(,I) FOR SPLINE
      CALL SPLNXX(RMINSV,DXSV,NPSV,X(IXY+1),X(IXSC+1))
      IQ=IQ+NQPL
      GO TO 9000
 
C     SKIP OVER POINTS IF WE ARE NOT RETAINING THIS SYMMETRY
 9001 READ(IIN,503,END=9100) (X(IXX+IX),IX=1,NP)
C
 9000 CONTINUE
C
      IF (NSKIP.GT.MXSKPR.OR.(NSKIP.GT.0.AND.IPRINT.LE.0))
     1       WRITE(6,604) NSKIP
  604 FORMAT('0 *** TOTAL OF',I6,' SYMMETRIES WERE SKIPPED ON &POTL',
     1       ' INPUT CRITERIA.')
C     FREE UP STORAGE USED FOR X(NPSV) AND NO LONGER NEEDED
      MX=MX+NPSV
C     SET PARAMETERS INTO CALLING ARGUMENTS
      NPOTL=MXLAM
C     BELOW FOR NEW (FEB 94) USE OF IV WITH BOTH ITYPE=2 AND 7
      IF (ITYPE.EQ.2.OR.ITYPE.EQ.7) THEN
        NPOTL=0
        DO 2010 I=1,MXLAM
 2010   NPOTL=MAX0(NPOTL,LAM((I-1)*NQPL+1))
        NPOTL=NPOTL+1
      ENDIF
C
      WRITE(6,661) EPSIL,RM,MXLAM,NPOTL
  661 FORMAT('0 ENERGY IN UNITS OF EPSILON =',F15.5,' CM-1'/
     1   '  R      IN UNITS OF RM      =',F15.5,' ANGSTROMS'/
     2   '0 MXLAM =',I5/'  NPOTL =',I5)
C
      R=RM
      P(1)=EPSIL
      MPOTL=NPOTL
      MXLMB=MXLAM
      RETURN
C
C     UNEXPECTED END OF FILE ON (IIN)
 9100 WRITE(6,*) ' *** POTENL(SPLINE).  UNEXPECTED EOF ON INPUT TAPE.'
      WRITE(6,*) ' ***                  TERMINATING.'
      STOP
      END
      SUBROUTINE SPLNXX (XMIN,DR,N,Y,SCR)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(N), SCR(N)
      DATA MAXTRY/100/
C
C     DO INITIALIZATION FOR CUBIC SPLINE FIT.
      NM1=N-1
      DR2=DR*DR
      H2=2D0*DR2
      OMEGA=4D0*(2D0-DSQRT(3D0))
      EPS=1D-8
C  DEFINITION OF 'NATURAL' SPLINE HAS S2(1) = S2(N) = 0.
C     NOT USED HERE.  ASSUME THAT SCR(1), SCR(N) SET IN CALLING PROG.
C     SCR(1)=0D0
C     SCR(N)=0D0
C
      DO 52 I=2,NM1
   52 SCR(I)=(Y(I+1)-2D0*Y(I)+Y(I-1))/DR2
      DO 521 ICNT=1,MAXTRY
      ETA=0.
      DO 10 I=2,NM1
       DELSQY=(Y(I+1)-2D0*Y(I)+Y(I-1))/H2
      W=(3D0*DELSQY-.25D0*(SCR(I-1)+SCR(I+1))-SCR(I))*OMEGA
      SCR(I)=SCR(I)+W
      IF (SCR(I).NE.0D0)  W=W/SCR(I)
   10 ETA=DMAX1(ETA,DABS(W))
C
      IF (ETA.GE.EPS)  GO TO 521
C     SCR(), WHICH IS EQUAL TO 2ND DERIV. AT ABSCISSAE, HAS CONVERGED.
C     WRITE(6,600)  ICNT,(SCR(I),I=1,N)
  600 FORMAT('0 SPLINE INIT. CONVERGED IN',I4,' ITERATIONS.'/
     1    (' ',8D16.8))
      RETURN
  521 CONTINUE
      WRITE(6,601)  MAXTRY
  601 FORMAT('0 * * * ERROR.  SPLINE INIT. DID NOT CONVERGE IN',I4
     1   ,' ITERATIONS FOR 2ND DERIVATIVE AT KNOTS.')
      STOP
      END
Modified: Wed Mar 8 17:00:00 1995 GMT
Page accessed 9179 times since Sat Apr 17 21:25:25 1999 GMT