molscat
|
README,
README.jkl,
README.v12,
dblas.f,
dblas.f.Z,
dcs_save.f,
diag_eispack.f,
ghm_save.f,
ghm_subs.f,
ghm_vib.f,
ident2disting.f,
lapack.f,
lapack.f.Z,
prbr_save.f,
prbr_vib.f,
read_isigu.f,
restrt.v12.f,
sbe.doc,
sbe_save.f,
sig_save.f,
spline.f,
syminv.f,
test1.input,
test1.v12.out,
test1.v14.out,
test2.input,
test2.v12.out,
test2.v14.out,
test3.f,
test3.input,
test3.v12.input,
test3.v12.out,
test3.v14.out,
test5.f,
test5.input,
test5.v12.out,
test5.v14.out,
test6.input,
test6.v12.out,
test6.v14.out,
test8.input,
test8.v12.out,
test8.v14.out,
timers.f,
timers_c.c,
v12.f,
v14.doc.tar,
v14.f,
v14.f.Z,
version_12.doc,
version_14.doc,
version_14.tutorial,
|
|
|
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
|