C C ******************************** C ** CENTRAL FIELD ORBITALS ** C ******************************** C C CALCULATION OF STATIONARY STATES OF ATOMIC CENTRAL FIELDS. C THE ADOPTED POTENTIALS COMBINE THE COULOMB FIELD OF A POINT OR C FINITE NUCLEUS WITH THE ANALYTICAL DHFS SCREENING FUNCTION OF C NEUTRAL ATOMS AND A LOCAL EXCHANGE CORRECTION (WHICH REMOVES C THE SELF-INTERACTION TERM FROM THE SCREENED FIELD). C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,NPPG=NDIM+1,NPTG=NDIM+NPPG) PARAMETER(HREV=27.211396D0) PARAMETER (PI=3.1415926535897932D0, PI4=4.0D0*PI) DIMENSION R(NDIM),RV(NDIM) C **** OUTPUT RADIAL FUNCTIONS. COMMON/RADWF/RAD(NDIM),P(NDIM),Q(NDIM),NGP,ILAST,IER C C **** READ POTENTIAL DATA. C 10 CONTINUE WRITE(6,*) ' ENTER ATOMIC NUMBER ...' READ(5,*) IZ IF(IZ.GT.100) STOP WRITE(6,*) ' ENTER MASS NUMBER ...' READ(5,*) ATW RNUC=2.3D-5*ATW**0.33333333333D0 WRITE(6,*) ' DHFS SCREENING? (1=YES, 0=NO) ...' READ(5,*) ISCR IF(ISCR.NE.0) ISCR=1 IF(ISCR.EQ.1) THEN WRITE(6,*) ' INCLUDE EXCHANGE POTENTIAL? (1=YES, 0=NO)', 1 ' ...' READ(5,*) ISINT IF(ISINT.NE.0) ISINT=1 ENDIF Z=IZ WRITE(6,'(3X,''POTENTIAL: '')') WRITE(6,'(3X,''ATOMIC NUMBER = '',I3)') IZ WRITE(6,'(3X,''MASS NUMBER = '',F7.3)') ATW WRITE(6,'(3X,''NUCLEAR RADIUS ='',1P,E12.5)') RNUC IF(ISCR.EQ.1) THEN WRITE(6,'(3X,''DHFS SCREENING'')') IF(ISINT.EQ.0) THEN WRITE(6,'(3X,''WITHOUT EXCHANGE POTENTIAL'')') ELSE WRITE(6,'(3X,''WITH EXCHANGE POTENTIAL''/)') ENDIF ELSE WRITE(6,'(3X,''BARE NUCLEUS, HIDROGENIC ION''/)') ENDIF C C C **** POTENTIAL GRID. C NV=NDIM NGP=NDIM CALL VATOM(IZ,RNUC,R,RV,NV,ISCR,ISINT) C OPEN(7,FILE='POTEN.DAT') RAD(1)=0.0D0 DO 11 I=2,NV RAD(I)=R(I) WRITE(7,'(1X,1P,3E15.7)') R(I),RV(I),RV(I)/R(I) 11 CONTINUE CLOSE(UNIT=7) C CALL VINT(R,RV,NV) C 20 CONTINUE WRITE(6,*) ' ' WRITE(6,*) ' SELECT ONE OPTION ...' WRITE(6,*) ' 1: SCHRODINGER EQUATION.' WRITE(6,*) ' 2: DIRAC EQUATION.' READ(5,*) IOPT WRITE(6,*) ' ' C C **** SCHRODINGER EQUATION. BOUND STATE. C IF(IOPT.EQ.1) THEN WRITE(6,*) ' ENTER N, L AND EPS ...' READ(5,*) N,L,EPS EPS=DMAX1(EPS,1.0D-15) IF(N.LT.1.OR.L.LT.0.OR.L.GE.N) STOP E=-Z**2/(2.0D0*N*N) CALL SBOUND(E,EPS,EPS,N,L) IF(IER.NE.0) GO TO 20 EEV=E*HREV WRITE(6,1100) N,L,EPS,E,EEV 1100 FORMAT(/1X,1P,'**** SCHRODINGER EQUATION. ', 2 /7X,'STATE: N=',I4,', L=',I4,' (EPS=',E8.1,')' 3 /7X,'ENERGY=',E22.14,' A.U.=',E22.14,' EV') C OPEN(7,FILE='SCHROD.DAT') WRITE(7,'(3X,''POTENTIAL: '')') WRITE(7,'(3X,''ATOMIC NUMBER = '',I3)') IZ WRITE(7,'(3X,''MASS NUMBER = '',F7.3)') ATW WRITE(7,'(3X,''NUCLEAR RADIUS ='',1P,E12.5)') RNUC IF(ISCR.EQ.1) THEN WRITE(7,'(3X,''DHFS SCREENING'')') IF(ISINT.EQ.0) THEN WRITE(7,'(3X,''WITHOUT EXCHANGE POTENTIAL'')') ELSE WRITE(7,'(3X,''WITH EXCHANGE POTENTIAL''/)') ENDIF ELSE WRITE(7,'(3X,''BARE NUCLEUS''/)') ENDIF WRITE(7,'(3X,''SCHRODINGER ENERGY LEVEL:'')') WRITE(7,'(3X,''N ='',I3)') N WRITE(7,'(3X,''L ='',I3)') L WRITE(7,'(3X,''EPS = '',1P,E9.2)') EPS WRITE(7,'(3X,''ENERGY = '',1P,E22.14,'' AU'')') E WRITE(7,'(3X,'' = '',1P,E22.14,'' eV''//)') EEV WRITE(7,'(8X,''R'',11X,''RDEN'',10X,''P'')') DO 30 I=2,NGP IF(DABS(P(I)).LT.1.0D-35) P(I)=1.0D-35 DEN=P(I)*P(I) WRITE(7,'(1X,1P,4E13.5)') RAD(I),DEN,P(I) 30 CONTINUE CLOSE(UNIT=7) OPEN(7,FILE='POTEN.DAT') WRITE(7,'(8X,''R'',14X,''V(R)'',8X,''V RADIAL'')') DO 31 I=2,NGP VEFF=RV(I)/R(I)+(L*(L+1.0D0))/(2.0D0*R(I)**2) WRITE(7,'(1X,1P,3E15.7)') R(I),RV(I),VEFF 31 CONTINUE CLOSE(UNIT=7) C C **** DIRAC EQUATION. BOUND STATE. C ELSE WRITE(6,*) ' ENTER N, L, J AND EPS ...' READ(5,*) N,L,RJ,EPS I=RJ+RJ+0.0001 RJ=0.5D0*I IF(DABS(L-RJ).GT.0.50000001D0) STOP EPS=DMAX1(EPS,1.0D-15) K=(L-RJ)*(RJ+RJ+1.00000001D0) IF(N.LT.1.OR.K.EQ.0.OR.K.GE.N.OR.K.LT.-N) STOP E=-Z**2/(2.0D0*N*N) CALL DBOUND(E,EPS,EPS,N,K) IF(IER.NE.0) GO TO 20 EEV=E*HREV WRITE(6,1200) N,L,RJ,EPS,E,EEV 1200 FORMAT(/1X,'**** DIRAC EQUATION. ',/7X,'STATE: N=',I4, 1 ', L=',I4,', J=',F6.1,' (EPS=',1PE8.1,')', 3 /7X,'ENERGY=',E22.14,' A.U.=',E22.14,' EV') C OPEN(7,FILE='DIRAC.DAT') WRITE(7,'(3X,''POTENTIAL: '')') WRITE(7,'(3X,''ATOMIC NUMBER = '',I3)') IZ WRITE(7,'(3X,''MASS NUMBER = '',F7.3)') ATW WRITE(7,'(3X,''NUCLEAR RADIUS ='',1P,E12.5)') RNUC IF(ISCR.EQ.1) THEN WRITE(7,'(3X,''DHFS SCREENING'')') IF(ISINT.EQ.0) THEN WRITE(7,'(3X,''WITHOUT EXCHANGE POTENTIAL'')') ELSE WRITE(7,'(3X,''WITH EXCHANGE POTENTIAL''/)') ENDIF ELSE WRITE(7,'(3X,''BARE NUCLEUS''/)') ENDIF WRITE(7,'(3X,''DIRAC ENERGY LEVEL:'')') WRITE(7,'(3X,''N ='',I3)') N WRITE(7,'(3X,''L ='',I3)') L WRITE(7,'(3X,''J ='',F5.1)') RJ WRITE(7,'(3X,''EPS = '',1P,E9.2)') EPS WRITE(7,'(3X,''ENERGY = '',1P,E22.14,'' AU'')') E WRITE(7,'(3X,'' = '',1P,E22.14,'' eV''//)') EEV WRITE(7,'(8X,''R'',11X,''RDEN'',10X,''P'',11X,''Q'')') DO 40 I=2,NGP IF(DABS(P(I)).LT.1.0D-35) P(I)=1.0D-35 IF(DABS(Q(I)).LT.1.0D-35) Q(I)=1.0D-35 DEN=P(I)*P(I)+Q(I)*Q(I) WRITE(7,'(1X,1P,4E13.5)') RAD(I),DEN,P(I),Q(I) 40 CONTINUE CLOSE(UNIT=7) OPEN(7,FILE='POTEN.DAT') WRITE(7,'(8X,''R'',14X,''V(R)'',8X,''V RADIAL'')') DO 41 I=2,NGP VEFF=RV(I)/R(I)+(L*(L+1.0D0))/(2.0D0*R(I)**2) WRITE(7,'(1X,1P,3E15.7)') R(I),RV(I),VEFF 41 CONTINUE CLOSE(UNIT=7) ENDIF GO TO 20 END C ************************************************************** C SUBROUTINE VATOM C ************************************************************** SUBROUTINE VATOM(IZ,RNUC,R,RV,NV,ISCR,ISINT) C C ANALYTICAL FIELD FOR THE CALCULATION OF ATOMIC ORBITALS. C THE POTENTIAL ENERGY CONSISTS OF 1) THE ELECTROSTATIC FIELD OF C A UNIFORMLY CHARGED SPHERICAL NUCLEUS OF RADIUS RNUC, 2) THE C ANALYTICAL DHFS SCREENING CORRECTION FOR NEUTRAL ATOMS AND C 3) A LOCAL EXCHANGE CORRECTION THAT ACOUNTS FOR THE SELF- C INTERACTION TERM IN THE ELECTRONIC INTERACTION POTENTIAL. C C INPUT PARAMETERS: C IZ ........ ATOMIC NUMBER. C RNUC ...... NUCLEAR RADIUS. C NV ........ NUMBER OF GRID POINTS (LARGER THAN 100). C ISCR ...... SCREENING CORRECTIONS ARE INCLUDED IN V(R) C IF ISCR=1. C ISINT ..... EXCHANGE CORRECTIONS ARE CONSIDERED WHEN C ISINT=1 AND ISCR=1. C C ** ALL QUANTITIES IN ATOMIC UNITS ** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,FR1O6=1.0D0/6.0D0) DIMENSION R(NDIM),RV(NDIM) IF(NV.LT.100) STOP C **** ANALYTICAL DHFS SCREENING FUNCTION. CALL DHFS(IZ,A1,A2,A3,AL1,AL2,AL3) Z=IZ C **** GRID POINTS. I=1 R(1)=0.0D0 IF(RNUC.GT.0.0D0) THEN INT=40 DR=RNUC/DFLOAT(INT) 1 I=I+1 R(I)=R(I-1)+DR IF(I.LE.INT+1) GO TO 1 ELSE I=2 R(2)=0.01D0/DMAX1(AL1,AL2,AL3) ENDIF FACT=DEXP(DLOG(2000.0D0/R(I))/(NV-I)) DO 2 J=I+1,NV R(J)=R(J-1)*FACT 2 CONTINUE C C **************** POTENTIAL ENERGY AT THE GRID POINTS. C ALPHA=1.5D0 IF(IZ.LT.3) ALPHA=0.0D0 FEX=(4.235654288D-1*ALPHA)**3 DO 3 I=1,NV RR=R(I) C **** NUCLEAR FIELD. IF(RR.LT.RNUC) THEN RVNUC=-0.5D0*Z*(3.0D0-(RR/RNUC)**2)*(RR/RNUC) ELSE RVNUC=-Z ENDIF C **** DHFS SCREENING FUNCTION. PHI=A1*DEXP(-AL1*RR)+A2*DEXP(-AL2*RR)+A3*DEXP(-AL3*RR) RVD=RVNUC*(PHI-1.0D0) C **** ELECTRON DENSITY AND LOCAL EXCHANGE FIELD. RHO=Z*(A1*AL1*AL1*DEXP(-AL1*RR)+A2*AL2*AL2*DEXP(-AL2*RR) 1 +A3*AL3*AL3*DEXP(-AL3*RR)) RVEX=-((FEX*RHO*RR*RR)**2+(RVD/Z)**6)**FR1O6 C RV(I)=RVNUC IF(ISCR.EQ.1) THEN RV(I)=RV(I)+RVD IF(ISINT.EQ.1) RV(I)=RV(I)+RVEX ENDIF 3 CONTINUE RETURN END C ************************************************************** C SUBROUTINE DHFS C ************************************************************** SUBROUTINE DHFS(IZ,A1,A2,A3,AL1,AL2,AL3) C C DHFS ANALYTICAL SCREENING FUNCTION PARAMETERS FOR FREE C NEUTRAL ATOMS. THE INPUT ARGUMENT IS THE ATOMIC NUMBER (IZ, C AN INTEGER, LESS THAN 104) C C REF.: F. SALVAT ET AL., PHYS. REV. A36 (1987) 467-474. C ELEMENTS FROM Z=93 TO 103 ADDED IN MARCH 1992. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION B1A(50),B1B(53),B2A(50),B2B(53),BL1A(50), 1 BL1B(53),BL2A(50),BL2B(53),BL3A(50),BL3B(53) DATA B1A /-7.05665D-6,-2.25920D-1,6.04537D-1,3.27766D-1, 1 2.32684D-1,1.53676D-1,9.95750D-2,6.25130D-2,3.68040D-2, 2 1.88410D-2,7.44440D-1,6.42349D-1,6.00152D-1,5.15971D-1, 3 4.38675D-1,5.45871D-1,7.24889D-1,2.19124D+0,4.85607D-2, 4 5.80017D-1,5.54340D-1,1.11950D-2,3.18350D-2,1.07503D-1, 5 4.97556D-2,5.11841D-2,5.00039D-2,4.73509D-2,7.70967D-2, 6 4.00041D-2,1.08344D-1,6.09767D-2,2.11561D-2,4.83575D-1, 7 4.50364D-1,4.19036D-1,1.73438D-1,3.35694D-2,6.88939D-2, 8 1.17552D-1,2.55689D-1,2.69313D-1,2.20138D-1,2.75057D-1, 9 2.71053D-1,2.78363D-1,2.56210D-1,2.27100D-1,2.49215D-1, A 2.15313D-1/ DATA B1B /1.80560D-1,1.30772D-1,5.88293D-2,4.45145D-1, B 2.70796D-1,1.72814D-1,1.94726D-1,1.91338D-1,1.86776D-1, C 1.66461D-1,1.62350D-1,1.58016D-1,1.53759D-1,1.58729D-1, D 1.45327D-1,1.41260D-1,1.37360D-1,1.33614D-1,1.29853D-1, E 1.26659D-1,1.28806D-1,1.30256D-1,1.38420D-1,1.50030D-1, F 1.60803D-1,1.72164D-1,1.83411D-1,2.23043D-1,2.28909D-1, G 2.09753D-1,2.70821D-1,2.37958D-1,2.28771D-1,1.94059D-1, H 1.49995D-1,9.55262D-2,3.19155D-1,2.40406D-1,2.26579D-1, I 2.17619D-1,2.41294D-1,2.44758D-1,2.46231D-1,2.55572D-1, J 2.53567D-1,2.43832D-1,2.41898D-1,2.44050D-1,2.40237D-1, K 2.34997D-1,2.32114D-1,2.27937D-1,2.29571D-1/ DATA B2A /-1.84386D+2,1.22592D+0,3.95463D-1,6.72234D-1, 1 7.67316D-1,8.46324D-1,9.00425D-1,9.37487D-1,9.63196D-1, 2 9.81159D-1,2.55560D-1,3.57651D-1,3.99848D-1,4.84029D-1, 3 5.61325D-1,-5.33329D-1,-7.54809D-1,-2.2852D0,7.75935D-1, 4 4.19983D-1,4.45660D-1,6.83176D-1,6.75303D-1,7.16172D-1, 5 6.86632D-1,6.99533D-1,7.14201D-1,7.29404D-1,7.95083D-1, 6 7.59034D-1,7.48941D-1,7.15671D-1,6.70932D-1,5.16425D-1, 7 5.49636D-1,5.80964D-1,7.25336D-1,7.81581D-1,7.20203D-1, 8 6.58088D-1,5.82051D-1,5.75262D-1,5.61797D-1,5.94338D-1, 9 6.11921D-1,6.06653D-1,6.50520D-1,6.15496D-1,6.43990D-1, A 6.11497D-1/ DATA B2B /5.76688D-1,5.50366D-1,5.48174D-1,5.54855D-1, B 6.52415D-1,6.84485D-1,6.38429D-1,6.46684D-1,6.55810D-1, C 7.05677D-1,7.13311D-1,7.20978D-1,7.28385D-1,7.02414D-1, D 7.42619D-1,7.49352D-1,7.55797D-1,7.61947D-1,7.68005D-1, E 7.73365D-1,7.52781D-1,7.32428D-1,7.09596D-1,6.87141D-1, F 6.65932D-1,6.46849D-1,6.30598D-1,6.17575D-1,6.11402D-1, G 6.00426D-1,6.42829D-1,6.30789D-1,6.21959D-1,6.10455D-1, H 6.03147D-1,6.05994D-1,6.23324D-1,6.56665D-1,6.42246D-1, I 6.24013D-1,6.30394D-1,6.29816D-1,6.31596D-1,6.49005D-1, J 6.53604D-1,6.43738D-1,6.48850D-1,6.70318D-1,6.76319D-1, K 6.65571D-1,6.88406D-1,6.94394D-1,6.82014D-1/ DATA BL1A /4.92969D+0,5.52725D+0,2.81741D+0,4.54302D+0, 1 5.99006D+0,8.04043D+0,1.08122D+1,1.48233D+1,2.14001D+1, 2 3.49994D+1,4.12050D+0,4.72663D+0,5.14051D+0,5.84918D+0, 3 6.67070D+0,6.37029D+0,6.21183D+0,5.54701D+0,3.02597D+1, 4 6.32184D+0,6.63280D+0,9.97569D+1,4.25330D+1,1.89587D+1, 5 3.18642D+1,3.18251D+1,3.29153D+1,3.47580D+1,2.53264D+1, 6 4.03429D+1,2.01922D+1,2.91996D+1,6.24873D+1,8.78242D+0, 7 9.33480D+0,9.91420D+0,1.71659D+1,5.52077D+1,3.13659D+1, 8 2.20537D+1,1.42403D+1,1.40442D+1,1.59176D+1,1.43137D+1, 9 1.46537D+1,1.46455D+1,1.55878D+1,1.69141D+1,1.61552D+1, A 1.77931D+1/ DATA BL1B /1.98751D+1,2.41540D+1,3.99955D+1,1.18053D+1, B 1.65915D+1,2.23966D+1,2.07637D+1,2.12350D+1,2.18033D+1, C 2.39492D+1,2.45984D+1,2.52966D+1,2.60169D+1,2.54973D+1, D 2.75466D+1,2.83460D+1,2.91604D+1,2.99904D+1,3.08345D+1, E 3.16806D+1,3.13526D+1,3.12166D+1,3.00767D+1,2.86302D+1, F 2.75684D+1,2.65861D+1,2.57339D+1,2.29939D+1,2.28644D+1, G 2.44080D+1,2.09409D+1,2.29872D+1,2.37917D+1,2.66951D+1, H 3.18397D+1,4.34890D+1,2.00150D+1,2.45012D+1,2.56843D+1, I 2.65542D+1,2.51930D+1,2.52522D+1,2.54271D+1,2.51526D+1, J 2.55959D+1,2.65567D+1,2.70360D+1,2.72673D+1,2.79152D+1, K 2.86446D+1,2.93353D+1,3.01040D+1,3.02650D+1/ DATA BL2A /2.00272D+0,2.39924D+0,6.62463D-1,9.85154D-1, 1 1.21347D+0,1.49129D+0,1.76868D+0,2.04035D+0,2.30601D+0, 2 2.56621D+0,8.71798D-1,1.00247D+0,1.01529D+0,1.17314D+0, 3 1.34102D+0,2.55169D+0,3.38827D+0,4.56873D+0,3.12426D+0, 4 1.00935D+0,1.10227D+0,4.12865D+0,3.94043D+0,3.06375D+0, 5 3.78110D+0,3.77161D+0,3.79085D+0,3.82989D+0,3.39276D+0, 6 3.94645D+0,3.47325D+0,4.12525D+0,4.95015D+0,1.69671D+0, 7 1.79002D+0,1.88354D+0,3.11025D+0,4.28418D+0,4.24121D+0, 8 4.03254D+0,2.97020D+0,2.86107D+0,3.36719D+0,2.73701D+0, 9 2.71828D+0,2.61549D+0,2.74124D+0,3.08408D+0,2.88189D+0, A 3.29372D+0/ DATA BL2B /3.80921D+0,4.61191D+0,5.91318D+0,1.79673D+0, B 2.69645D+0,3.45951D+0,3.46574D+0,3.48193D+0,3.50982D+0, C 3.51987D+0,3.55603D+0,3.59628D+0,3.63834D+0,3.73639D+0, D 3.72882D+0,3.77625D+0,3.82444D+0,3.87344D+0,3.92327D+0, E 3.97271D+0,4.09040D+0,4.20492D+0,4.24918D+0,4.24261D+0, F 4.23412D+0,4.19992D+0,4.14615D+0,3.73461D+0,3.69138D+0, G 3.96429D+0,3.24563D+0,3.62172D+0,3.77959D+0,4.25824D+0, H 4.92848D+0,5.85205D+0,2.90906D+0,3.55241D+0,3.79223D+0, I 4.00437D+0,3.67795D+0,3.63966D+0,3.61328D+0,3.43021D+0, J 3.43474D+0,3.59089D+0,3.59411D+0,3.48061D+0,3.50331D+0, K 3.61870D+0,3.55697D+0,3.58685D+0,3.64085D+0/ DATA BL3A /1.99732D+0,1.00000D+0,1.00000D+0,1.00000D+0, 1 1.00000D+0,1.00000D+0,1.00000D+0,1.00000D+0,1.00000D+0, 2 1.00000D+0,1.00000D+0,1.00000D+0,1.00000D+0,1.00000D+0, 3 1.00000D+0,1.67534D+0,1.85964D+0,2.04455D+0,7.32637D-1, 4 1.00000D+0,1.00000D+0,1.00896D+0,1.05333D+0,1.00137D+0, 5 1.12787D+0,1.16064D+0,1.19152D+0,1.22089D+0,1.14261D+0, 6 1.27594D+0,1.00643D+0,1.18447D+0,1.35819D+0,1.00000D+0, 7 1.00000D+0,1.00000D+0,7.17673D-1,8.57842D-1,9.47152D-1, 8 1.01806D+0,1.01699D+0,1.05906D+0,1.15477D+0,1.10923D+0, 9 1.12336D+0,1.43183D+0,1.14079D+0,1.26189D+0,9.94156D-1, A 1.14781D+0/ DATA BL3B /1.28288D+0,1.41954D+0,1.54707D+0,1.00000D+0, B 6.81361D-1,8.07311D-1,8.91057D-1,9.01112D-1,9.10636D-1, C 8.48620D-1,8.56929D-1,8.65025D-1,8.73083D-1,9.54998D-1, D 8.88981D-1,8.96917D-1,9.04803D-1,9.12768D-1,9.20306D-1, E 9.28838D-1,1.00717D+0,1.09456D+0,1.16966D+0,1.23403D+0, F 1.29699D+0,1.35350D+0,1.40374D+0,1.44284D+0,1.48856D+0, G 1.53432D+0,1.11214D+0,1.23735D+0,1.25338D+0,1.35772D+0, H 1.46828D+0,1.57359D+0,7.20714D-1,8.37599D-1,9.33468D-1, I 1.02385D+0,9.69895D-1,9.82474D-1,9.92527D-1,9.32751D-1, J 9.41671D-1,1.01827D+0,1.02554D+0,9.66447D-1,9.74347D-1, K 1.04137D+0,9.90568D-1,9.98878D-1,1.04473D+0/ C IIZ=IABS(IZ) IF(IIZ.GT.103) IIZ=103 IF(IIZ.EQ.0) IIZ=1 IF(IIZ.LE.50) THEN A1=B1A(IIZ) A2=B2A(IIZ) AL1=BL1A(IIZ) AL2=BL2A(IIZ) AL3=BL3A(IIZ) ELSE A1=B1B(IIZ-50) A2=B2B(IIZ-50) AL1=BL1B(IIZ-50) AL2=BL2B(IIZ-50) AL3=BL3B(IIZ-50) ENDIF A3=1.0D0-(A1+A2) RETURN END C C ***************************************** C *** SUBROUTINE PACKAGE RADIAL *** C ***************************************** C C BARCELONA, DECEMBER 1994. C C NUMERICAL SOLUTION OF THE SCHRODINGER (S) AND DIRAC (D) RADIAL C WAVE EQUATIONS. CUBIC SPLINE FIELD + POWER SERIES METHOD. C C IT IS ASSUMED THAT THE POTENTIAL ENERGY V(R) IS SUCH THAT THE C FUNCTION R*V(R) IS FINITE FOR ALL R AND TENDS TO CONSTANT VA- C LUES WHEN R TENDS TO ZERO AND TO INFINITY. C C**** ALL QUANTITIES ARE IN ATOMIC HARTREE UNITS. C FOR ELECTRONS AND POSITRONS C UNIT OF LENGTH = A0 = 5.29177D-11 METRES (= BOHR RADIUS), C UNIT OF ENERGY = E0 = 27.2114 EV (= HARTREE ENERGY). C FOR PARTICLES OF MASS 'M' (IN UNITS OF THE ELECTRON MASS) THE C CORRESPONDING UNITS ARE C UNIT OF LENGTH = A0/M, C UNIT OF ENERGY = M*E0. C C C THE CALLING SEQUENCE FROM THE MAIN PROGRAM IS: C C**** CALL VINT(R,RV,NV) C C THIS IS AN INITIALIZATION ROUTINE. IT DETERMINES THE NATURAL C CUBIC SPLINE THAT INTERPOLATES THE TABLE OF VALUES OF THE C FUNCTION R*V(R) PROVIDED BY THE USER. C INPUT ARGUMENTS: C R(I) ..... INPUT POTENTIAL GRID POINTS (REPEATED VALUES ARE C INTERPRETED AS DISCONTINUITIES). C RV(I) .... R(I) TIMES THE POTENTIAL ENERGY AT R=R(I). C NV ....... NUMBER OF POINTS IN THE TABLE (.LE.NDIM). C THE R(I) GRID _MUST_ INCLUDE THE ORIGIN (R=0), AND EXTEND UP C TO RADIAL DISTANCES FOR WHICH THE FUNCTION R*V(R) REACHES ITS C (CONSTANT) ASYMPTOTIC VALUE. THE INPUT GRID POINTS MUST BE IN C NON-DECREASING ORDER, I.E. R(I+1).GE.R(I). C C**** CALL SBOUND(E,EPS,DELL,N,L) OR DBOUND(E,EPS,DELL,N,K) C C THESE SUBROUTINES SOLVE THE RADIAL WAVE EQUATIONS FOR BOUND C STATES. C INPUT ARGUMENTS: C E ........ ESTIMATED BINDING ENERGY (A GOOD INITIAL ESTI- C MATE SPEEDS UP THE CALCULATION). C EPS ...... GLOBAL TOLERANCE, I.E. ALLOWED RELATIVE ERROR IN C THE SUMMATION OF THE RADIAL FUNCTION SERIES. C DELL ..... EIGENVALUE TOLERANCE, I.E. THE RELATIVE ERROR OF C THE COMPUTED EIGENVALUE IS REQUIRED TO BE LESS C THAN DELL (A CONVENIENT VALUE IS DELL=EPS). C N ........ PRINCIPAL QUANTUM NUMBER. C L ........ ORBITAL ANGULAR MOMENTUM QUANTUM NUMBER. C K ........ RELATIVISTIC ANGULAR MOMENTUM QUANTUM NUMBER. C (NOTE: 0.LE.L.LE.N-1, -N.LE.K.LE.N-1) C OUTPUT ARGUMENT: C E ........ BINDING ENERGY. C C**** CALL SFREE(E,EPS,PHASE,L) OR DFREE(E,EPS,PHASE,K) C C THESE SUBROUTINES SOLVE THE RADIAL WAVE EQUATIONS FOR FREE C STATES. C INPUT ARGUMENTS: C E ........ KINETIC ENERGY. C EPS ...... GLOBAL TOLERANCE, I.E. ALLOWED RELATIVE ERROR IN C THE SUMMATION OF THE RADIAL FUNCTION SERIES. C L ........ ORBITAL ANGULAR MOMENTUM QUANTUM NUMBER. C K ........ RELATIVISTIC ANGULAR MOMENTUM QUANTUM NUMBER. C (NOTE: L.GE.0, K.NE.0) C OUTPUT ARGUMENTS: C PHASE .... INNER PHASE SHIFT (IN RADIANS), CAUSED BY THE C SHORT RANGE COMPONENT OF THE FIELD. FOR MODIFIED C COULOMB FIELDS, WE HAVE C TOTAL PHASE SHIFT = PHASE + DELTA C WHERE DELTA IS THE COULOMB PHASE SHIFT (DELIVER- C ED THROUGH THE COMMON BLOCK C COMMON/OCOUL/RK,ETA,DELTA). C C C THE VALUES OF THE RADIAL FUNCTIONS ARE DELIVERED THROUGH THE C COMMON BLOCK C COMMON/RADWF/RAD(NDIM),P(NDIM),Q(NDIM),NGP,ILAST,IER C WITH NDIM=800 (CHANGE THE VALUE OF THIS PARAMETER IN _ALL_ C SUBROUTINES IF A LARGER NUMBER OF GRID POINTS IS NEEDED). THE C GRID OF POINTS RAD(I), WHERE THE RADIAL WAVE FUNCTIONS ARE C TABULATED, CAN BE ARBITRARILY SELECTED BY THE USER. THE QUAN- C TITIES C RAD(I) ... USER'S RADIAL GRID, C NGP ...... NUMBER OF GRID POINTS (.LE.NDIM), C MUST BE DEFINED BEFORE CALLING THE SOLUTION SUBROUTINES. THE C RADIAL POINTS ARE SORTED IN INCREASING ORDER BY THE SOLUTION C SUBROUTINES (TO AVOID CONFUSION, IT IS ADVISABLE TO ORDER THE C INPUT GRID IN THIS WAY). C THE OUTPUT QUANTITIES ARE C P(I) ..... VALUE OF THE RADIAL FUNCTION P(R) AT THE I-TH C GRID POINT, R=RAD(I). C Q(I) ..... VALUE OF THE RADIAL FUNCTION Q(R) AT THE I-TH C GRID POINT (= P'(R) FOR SCHRODINGER PARTICLES). C ILAST .... *** BOUND STATES: FOR R.GT.RAD(ILAST), P(R) AND C Q(R) ARE SET EQUAL TO 0.0D0. C *** FREE STATES: FOR R.GT.RAD(ILAST), P(R) AND C Q(R) ARE OBTAINED IN TERMS OF THE REGULAR AND C IRREGULAR ASYMPTOTIC COULOMB FUNCTIONS AS C P(R)=DCOS(PHASE)*FU(R)+DSIN(PHASE)*GU(R) C Q(R)=DCOS(PHASE)*FL(R)+DSIN(PHASE)*GL(R) C WHERE FU, GU AND FL, GL ARE GIVEN BY SUBROUTINES C SCOUL AND DCOUL WITH Z=RV(NV). WHEN THE ABSOLUTE C VALUE OF RV(NV) IS LESS THAN EPS, Z IS SET EQUAL C TO ZERO, SO THAT THE FUNCTIONS FU, GU, FL AND GL C REDUCE TO SPHERICAL BESSEL FUNCTIONS OF INTEGER C ORDER (GIVEN BY FUNCTION BESJN). C IER ...... ERROR CODE. A VALUE LARGER THAN ZERO INDICATES C THAT SOME FATAL ERROR HAS BEEN FOUND DURING THE C CALCULATION. C C C BOUND STATE WAVE FUNCTIONS ARE NORMALIZED TO UNITY. THE ADOPT- C ED NORMALIZATION FOR FREE STATES IS SUCH THAT P(R) OSCILLATES C WITH UNIT AMPLITUDE IN THE ASYMPTOTIC REGION. C C C**** ERROR CODES (AND TENTATIVE SOLUTIONS...): C 0 .... EVERYTHING IS O.K. C 1 .... EMIN.GE.0 IN 'BOUND' (USE A DENSER GRID. IF THE ERROR C PERSISTS THEN PROBABLY SUCH A BOUND STATE DOES NOT C EXIST). C 2 .... E=0 IN 'BOUND' (PROBABLY THIS BOUND STATE DOES NOT C EXIST). C 3 .... RAD(NGP) TOO SMALL IN 'BOUND' (EXTEND THE GRID TO C LARGER RADII). C 4 .... SEVERAL ZEROS OF P(R) IN A SINGLE INTERVAL IN 'BOUND' C (USE A DENSER GRID). C 5 .... E OUT OF RANGE IN 'BOUND' (ACCUMULATED ROUND-OFF C ERRORS?). C 6 .... RV(NGP)<<0 OR E>0 IN 'BOUND' (CHECK THE INPUT POTEN- C TIAL VALUES). C 7 .... E.LE.0 IN 'FREE'. C 8 .... RAD(NGP) TOO SMALL IN 'FREE' (EXTEND THE GRID TO C LARGER RADII). C C THE PROGRAM STOPS WHEN THE INPUT QUANTUM NUMBERS ARE OUT OF C RANGE. C C C**** YOU CAN GET A PRINTED MANUAL OF THE SUBROUTINE PACKAGE BY C RUNNING THE INPUT FILE 'RADIAL.TEX' THROUGH LaTeX, AND PRINT- C ING THE OUTPUT FILE. BRACKETED NUMBERS IN COMMENT LINES OF THE C PRESENT FILE INDICATE THE CORRESPONDING EQUATIONS IN THAT MA- C NUAL. C C ************************************************************** C SUBROUTINE VINT C ************************************************************** SUBROUTINE VINT(R,RV,NV) C C NATURAL CUBIC SPLINE INTERPOLATION FOR R*V(R) FROM THE C INPUT RADII AND POTENTIAL VALUES (128). C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,NPPG=NDIM+1,NPTG=NDIM+NPPG) COMMON/VGRID/RG(NPPG),RVG(NPPG),VA(NPPG),VB(NPPG),VC(NPPG), 1 VD(NPPG),NVT COMMON/RGRID/X(NPTG),P(NPTG),Q(NPTG),IND(NPTG),NRT COMMON/STORE/Y(NPTG),A(NPTG),B(NPTG),C(NPTG),D(NPTG) DIMENSION R(NDIM),RV(NDIM) C IF(R(1).LT.0.0D0) THEN WRITE(6,2101) 2101 FORMAT(1X,'*** ERROR IN VINT: R(1).LT.0.') STOP ENDIF IF(NV.GT.NDIM) THEN WRITE(6,2102) NDIM 2102 FORMAT(1X,'*** ERROR IN VINT: INPUT POTENTIAL GRID WITH ', 1 'MORE THAN ',I5,' DATA POINTS.') STOP ENDIF R(1)=0.0D0 C IO=0 I=0 K=0 1 I=I+1 K=K+1 X(K)=R(I) Y(K)=RV(I) IF(I.EQ.NV) GO TO 2 IF(R(I).LT.R(I+1)-1.0D-12) GO TO 1 2 CONTINUE C CALL SPLINE(X,Y,A,B,C,D,0.0D0,0.0D0,K) C K=K-1 DO 3 J=1,K IO=IO+1 RG(IO)=X(J) RVG(IO)=Y(J) VA(IO)=A(J) VB(IO)=B(J) VC(IO)=C(J) VD(IO)=D(J) 3 CONTINUE IF(I.LT.NV) THEN K=0 GO TO 1 ENDIF C **** AN EXTRA POINT IS ADDED TO THE GRID, AND R*V(R) IS SET C EQUAL TO RV(NV) FOR R.GE.R(NV) IO=IO+1 RG(IO)=X(K+1) RVG(IO)=Y(K+1) VA(IO)=RVG(IO) VB(IO)=0.0D0 VC(IO)=0.0D0 VD(IO)=0.0D0 NVT=IO+1 RG(NVT)=2.0D0*RG(IO) RVG(NVT)=RVG(IO) VA(NVT)=RVG(IO) VB(NVT)=0.0D0 VC(NVT)=0.0D0 VD(NVT)=0.0D0 RETURN END C ************************************************************** C SUBROUTINE SBOUND C ************************************************************** SUBROUTINE SBOUND(E,EPS,DELL,N,L) C C THIS SUBROUTINE SOLVES THE SCHRODINGER EQUATION FOR BOUND C STATES. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,NPPG=NDIM+1,NPTG=NDIM+NPPG) C **** SET IWR=0 TO AVOID PRINTING PARTIAL RESULTS. PARAMETER (IWR=1) COMMON/RADWF/RAD(NDIM),PIO(NDIM),QIO(NDIM),NGP,ILAST,IER COMMON/RGRID/R(NPTG),P(NPTG),Q(NPTG),IND(NPTG),NRT COMMON/VGRID/RG(NPPG),RV(NPPG),VA(NPPG),VB(NPPG),VC(NPPG), 1 VD(NPPG),NVT COMMON/STORE/Y(NPTG),A(NPTG),B(NPTG),C(NPTG),D(NPTG) COMMON/NZT/NZMAX C IER=0 IF(N.LT.1) THEN WRITE(6,2101) 2101 FORMAT(1X,'*** ERROR IN SBOUND: N.LT.1.') STOP ENDIF IF(L.LT.0) THEN WRITE(6,2102) 2102 FORMAT(1X,'*** ERROR IN SBOUND: L.LT.0.') STOP ENDIF IF(L.GE.N) THEN WRITE(6,2103) 2103 FORMAT(1X,'*** ERROR IN SBOUND: L.GE.N.') STOP ENDIF C **** RADIAL QUANTUM NUMBER. NR=N-L-1 C DELL=DABS(DELL) IF(E.GT.-1.0D-1) E=-1.0D-1 FL1=0.5D0*L*(L+1) C C **** MERGE THE 'RG' AND 'RAD' GRIDS, C T=DMAX1(0.5D0*EPS,1.0D-10) DO 100 I=1,NVT R(I)=RG(I) IND(I)=I 100 CONTINUE NRT=NVT DO 200 I=1,NGP RLOC=RAD(I) CALL FINDI(RG,RLOC,NVT,J) TST=DMIN1(DABS(RLOC-RG(J)),DABS(RLOC-RG(J+1))) IF(TST.LT.T) GO TO 200 NRT=NRT+1 R(NRT)=RLOC IND(NRT)=J 200 CONTINUE C **** ... AND SORT THE RESULTING R-GRID IN INCREASING ORDER. DO 400 I=1,NRT-1 RMIN=1.0D35 IMIN=I DO 300 J=I,NRT IF(R(J).LT.RMIN) THEN RMIN=R(J) IMIN=J ENDIF 300 CONTINUE IF(IMIN.NE.I) THEN RMIN=R(I) R(I)=R(IMIN) R(IMIN)=RMIN INDMIN=IND(I) IND(I)=IND(IMIN) IND(IMIN)=INDMIN ENDIF 400 CONTINUE C C **** MINIMUM OF THE EFFECTIVE RADIAL POTENTIAL. (165) C EMIN=1.0D0 DO 1 I=2,NRT RN=R(I) J=IND(I) RVN=VA(J)+RN*(VB(J)+RN*(VC(J)+RN*VD(J))) EMIN=DMIN1(EMIN,(RVN+FL1/RN)/RN) 1 CONTINUE IF(EMIN.GT.-1.0D-35) THEN IER=1 WRITE(6,1001) 1001 FORMAT(1X,'*** ERROR 1 (IN SBOUND): EMIN.GE.0.'/5X, 1'(USE A DENSER GRID. IF THE ERROR PERSISTS THEN PROBABLY', 2/6X,'SUCH A BOUND STATE DOES NOT EXIST).') RETURN ENDIF RN=R(NRT) J=IND(NRT) RVN=VA(J)+RN*(VB(J)+RN*(VC(J)+RN*VD(J))) ESUP=(RVN+FL1/RN)/RN IF(ESUP.GT.0.0D0) ESUP=0.0D0 IF(L.EQ.0) THEN TEHYDR=-(R(1)/N)**2 EMIN=DMIN1(E,TEHYDR,-10.0D0) ENDIF IF(E.GT.ESUP.OR.E.LT.EMIN) E=0.5D0*(ESUP+EMIN) EMAX=ESUP ICMIN=0 ICMAX=0 C 2 CONTINUE IF(E.GT.-1.0D-16) THEN IER=2 WRITE(6,1002) 1002 FORMAT(1X,'*** ERROR 2 (IN SBOUND): E=0.'/5X, 1'(PROBABLY THIS BOUND STATE DOES NOT EXIST).') RETURN ENDIF C **** OUTER TURNING POINT. (165) DO 3 K=2,NRT IOTP=NRT+2-K RN=R(IOTP) J=IND(IOTP) RVN=VA(J)+RN*(VB(J)+RN*(VC(J)+RN*VD(J))) EKIN=E-(RVN+FL1/RN)/RN IF(EKIN.GT.0.0D0) GO TO 4 3 CONTINUE 4 IOTP=IOTP+1 IF(IOTP.GT.NRT-1) THEN IER=3 WRITE(6,1003) 1003 FORMAT(1X,'*** ERROR 3 (IN SBOUND): RAD(NGP) TOO SMALL.' 1/5X,'(EXTEND THE GRID TO LARGER RADII).') RETURN ENDIF C C **** OUTWARD SOLUTION. C CALL SOUTW(E,EPS,L,NR,NZERO,IOTP) IF(NZMAX.GT.1.AND.NZERO.LE.NR) THEN IER=4 WRITE(6,1004) 1004 FORMAT(1X,'*** ERROR 4 (IN SBOUND): SEVERAL ZEROS OF P(R)', 1/5X,'IN A SINGLE INTERVAL (USE A DENSER GRID).') RETURN ENDIF C C **** TOO MANY NODES. IF(NZERO.GT.NR) THEN IF(ICMIN.EQ.0) EMIN=EMIN-2.0D0*(EMAX-EMIN) EMAX=E ICMAX=1 E=0.5D0*(EMIN+E) IF(IWR.NE.0) THEN WRITE(6,2000) N,L WRITE(6,2001) NR,NZERO,IOTP,NRT WRITE(6,2002) E WRITE(6,2004) EMIN,EMAX ENDIF C IF(EMAX-EMIN.LT.DELL*DABS(EMIN)) THEN IER=5 WRITE(6,1005) RETURN ENDIF C GO TO 2 ENDIF C **** TOO FEW NODES. IF(NZERO.LT.NR) THEN ICMIN=1 EMIN=E E=0.5D0*(E+EMAX) IF(IWR.NE.0) THEN WRITE(6,2000) N,L WRITE(6,2001) NR,NZERO,IOTP,NRT WRITE(6,2002) E WRITE(6,2004) EMIN,EMAX ENDIF C IF(EMAX-EMIN.LT.DELL*DABS(EMIN)) THEN IER=5 WRITE(6,1005) RETURN ENDIF C GO TO 2 ENDIF C **** THE CORRECT NUMBER OF NODES HAS BEEN FOUND. PO=P(IOTP) QO=Q(IOTP) C C **** INWARD SOLUTION C CALL SINW(E,EPS,L,IOTP) IF(IER.GT.0) RETURN C **** MATCHING OF THE OUTWARD AND INWARD SOLUTIONS. FACT=PO/P(IOTP) DO 5 I=IOTP,ILAST P(I)=P(I)*FACT 5 Q(I)=Q(I)*FACT QI=Q(IOTP) RLAST=R(ILAST) C **** NORMALIZATION CALL SPLINE(R,P,A,B,C,D,0.0D0,0.0D0,ILAST) CALL INTEG2(R,A,B,C,D,0.0D0,RLAST,SUM,ILAST) C **** EIGENVALUE CORRECTION. (175) IF(SUM.LT.1.0D-15) SUM=1.0D0 DE=PO*(QO-QI)/(SUM+SUM) EP=E+DE C IF(DE.LT.0.0D0) THEN ICMAX=1 EMAX=E ENDIF C IF(DE.GT.0.0D0) THEN ICMIN=1 EMIN=E ENDIF C IF(ICMIN.EQ.0.AND.EP.LT.EMIN) THEN EMIN=1.1D0*EMIN IF(EP.LT.EMIN) EP=0.5D0*(E+EMIN) ENDIF IF(ICMIN.EQ.1.AND.EP.LT.EMIN) EP=0.5D0*(E+EMIN) IF(ICMAX.EQ.1.AND.EP.GT.EMAX) EP=0.5D0*(E+EMAX) IF(EP.GT.ESUP) EP=0.5D0*(ESUP+E) C IF(IWR.NE.0) THEN WRITE(6,2000) N,L 2000 FORMAT(/2X,'SUBROUTINE SBOUND. N =',I3,' L =',I3) WRITE(6,2001) NR,NZERO,IOTP,NRT 2001 FORMAT(2X,'NR =',I3,' NZERO =',I3,' IOTP = ',I5, 1' NGP =',I5) WRITE(6,2002) EP 2002 FORMAT(2X,'E NEW = ',1P,D22.15) WRITE(6,2003) E,DE 2003 FORMAT(2X,'E OLD = ',1P,D22.15,' DE = ',D11.4) WRITE(6,2004) EMIN,EMAX 2004 FORMAT(2X,'EMIN = ',1P,D12.5,' EMAX = ',D12.5) ENDIF C IF(EP.GE.ESUP.AND.DABS(E-ESUP).LT.DELL*DABS(ESUP)) THEN IER=5 WRITE(6,1005) 1005 FORMAT(1X,'*** ERROR 5 (IN SBOUND): E OUT OF RANGE.'/5X, 1'(ACCUMULATED ROUND-OFF ERRORS?).') RETURN ENDIF EO=E E=EP IF(DMIN1(DABS(DE),DABS(E-EO)).GT.DABS(E*DELL)) GO TO 2 C **** NORMALIZATION FACT=1.0D0/DSQRT(SUM) DO 6 I=1,ILAST P(I)=P(I)*FACT 6 Q(I)=Q(I)*FACT C C **** EXTRACT THE 'RAD' GRID... C DO 500 I=1,NGP RLOC=RAD(I) CALL FINDI(R,RLOC,NRT,J) IF(RLOC-R(J).LT.R(J+1)-RLOC) THEN PIO(I)=P(J) QIO(I)=Q(J) ELSE PIO(I)=P(J+1) QIO(I)=Q(J+1) ENDIF 500 CONTINUE C IF(DABS(P(NRT)).GT.1.0D-5*DABS(P(IOTP))) THEN IER=3 WRITE(6,1003) ENDIF C RLOC=R(ILAST) CALL FINDI(RAD,RLOC,NGP,ILAST) ILAST=ILAST+1 RETURN END C ************************************************************** C SUBROUTINE DBOUND C ************************************************************** SUBROUTINE DBOUND(E,EPS,DELL,N,K) C C THIS SUBROUTINE SOLVES THE DIRAC EQUATION FOR BOUND STATES. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,NPPG=NDIM+1,NPTG=NDIM+NPPG, 1 SL=137.036D0) C **** SET IWR=0 TO AVOID PRINTING PARTIAL RESULTS. PARAMETER (IWR=1) COMMON/RADWF/RAD(NDIM),PIO(NDIM),QIO(NDIM),NGP,ILAST,IER COMMON/RGRID/R(NPTG),P(NPTG),Q(NPTG),IND(NPTG),NRT COMMON/VGRID/RG(NPPG),RV(NPPG),VA(NPPG),VB(NPPG),VC(NPPG), 1 VD(NPPG),NVT COMMON/STORE/Y(NPTG),A(NPTG),B(NPTG),C(NPTG),D(NPTG) COMMON/NZT/NZMAX C IER=0 IF(N.LT.1) THEN WRITE(6,2101) 2101 FORMAT(1X,'*** ERROR IN DBOUND: N.LT.1.') STOP ENDIF IF(K.EQ.0) THEN WRITE(6,2102) 2102 FORMAT(1X,'*** ERROR IN DBOUND: K.EQ.0.') STOP ENDIF IF(K.LT.-N) THEN WRITE(6,2103) 2103 FORMAT(1X,'*** ERROR IN DBOUND: K.LT.-N.') STOP ENDIF IF(K.GE.N) THEN WRITE(6,2104) 2104 FORMAT(1X,'*** ERROR IN DBOUND: K.GE.N.') STOP ENDIF C **** ORBITAL ANGULAR MOMENTUM QUANTUM NUMBER. (15) IF(K.LT.0) THEN L=-K-1 ELSE L=K ENDIF C **** RADIAL QUANTUM NUMBER. NR=N-L-1 C DELL=DABS(DELL) IF(E.GT.-1.0D-1) E=-1.0D-1 FL1=0.5D0*L*(L+1) C C **** MERGE THE 'RG' AND 'RAD' GRIDS, C T=DMAX1(0.5D0*EPS,1.0D-10) DO 100 I=1,NVT R(I)=RG(I) IND(I)=I 100 CONTINUE NRT=NVT DO 200 I=1,NGP RLOC=RAD(I) CALL FINDI(RG,RLOC,NVT,J) TST=DMIN1(DABS(RLOC-RG(J)),DABS(RLOC-RG(J+1))) IF(TST.LT.T) GO TO 200 NRT=NRT+1 R(NRT)=RLOC IND(NRT)=J 200 CONTINUE C **** ... AND SORT THE RESULTING R-GRID IN INCREASING ORDER. DO 400 I=1,NRT-1 RMIN=1.0D35 IMIN=I DO 300 J=I,NRT IF(R(J).LT.RMIN) THEN RMIN=R(J) IMIN=J ENDIF 300 CONTINUE IF(IMIN.NE.I) THEN RMIN=R(I) R(I)=R(IMIN) R(IMIN)=RMIN INDMIN=IND(I) IND(I)=IND(IMIN) IND(IMIN)=INDMIN ENDIF 400 CONTINUE C C **** MINIMUM OF THE EFFECTIVE RADIAL POTENTIAL. (165) C EMIN=1.0D0 DO 1 I=2,NRT RN=R(I) J=IND(I) RVN=VA(J)+RN*(VB(J)+RN*(VC(J)+RN*VD(J))) EMIN=DMIN1(EMIN,(RVN+FL1/RN)/RN) 1 CONTINUE IF(EMIN.GT.-1.0D-35) THEN IER=1 WRITE(6,1001) 1001 FORMAT(1X,'*** ERROR 1 (IN DBOUND): EMIN.GE.0.'/5X, 1'(USE A DENSER GRID. IF THE ERROR PERSISTS THEN PROBABLY', 2/6X,'SUCH A BOUND STATE DOES NOT EXIST).') RETURN ENDIF RN=R(NRT) J=IND(NRT) RVN=VA(J)+RN*(VB(J)+RN*(VC(J)+RN*VD(J))) ESUP=(RVN+FL1/RN)/RN IF(ESUP.GT.0.0D0) ESUP=0.0D0 IF(L.EQ.0) THEN TEHYDR=-(R(1)/N)**2 EMIN=DMIN1(E,TEHYDR,-10.0D0) ENDIF IF(EMIN.LT.-SL*SL) EMIN=-SL*SL IF(E.GT.ESUP.OR.E.LT.EMIN) E=0.5D0*(ESUP+EMIN) EMAX=ESUP ICMIN=0 ICMAX=0 C 2 CONTINUE IF(E.GT.-1.0D-16) THEN IER=2 WRITE(6,1002) 1002 FORMAT(1X,'*** ERROR 2 (IN DBOUND): E=0.'/5X, 1'(PROBABLY THIS BOUND STATE DOES NOT EXIST).') RETURN ENDIF C **** OUTER TURNING POINT. (165) DO 3 J=2,NRT IOTP=NRT+2-J RN=R(IOTP) JJ=IND(IOTP) RVN=VA(JJ)+RN*(VB(JJ)+RN*(VC(JJ)+RN*VD(JJ))) EKIN=E-(RVN+FL1/RN)/RN IF(EKIN.GT.0.0D0) GO TO 4 3 CONTINUE 4 IOTP=IOTP+1 IF(IOTP.GT.NRT-1) THEN IER=3 WRITE(6,1003) 1003 FORMAT(1X,'*** ERROR 3 (IN DBOUND): RAD(NGP) TOO SMALL.' 1/5X,'(EXTEND THE GRID TO LARGER RADII).') RETURN ENDIF C C **** OUTWARD SOLUTION. C CALL DOUTW(E,EPS,K,NR,NZERO,IOTP) IF(NZMAX.GT.1.AND.NZERO.LE.NR) THEN IER=4 WRITE(6,1004) 1004 FORMAT(1X,'*** ERROR 4 (IN DBOUND): SEVERAL ZEROS OF P(R)', 1/5X,'IN A SINGLE INTERVAL (USE A DENSER GRID).') RETURN ENDIF C C **** TOO MANY NODES. IF(NZERO.GT.NR) THEN IF(ICMIN.EQ.0) EMIN=EMIN-2.0D0*(EMAX-EMIN) EMAX=E ICMAX=1 E=0.5D0*(EMIN+E) IF(IWR.NE.0) THEN WRITE(6,2000) N,K WRITE(6,2001) NR,NZERO,IOTP,NRT WRITE(6,2002) E WRITE(6,2004) EMIN,EMAX ENDIF C IF(EMAX-EMIN.LT.DELL*DABS(EMIN)) THEN IER=5 WRITE(6,1005) RETURN ENDIF C GO TO 2 ENDIF C **** TOO FEW NODES. IF(NZERO.LT.NR) THEN ICMIN=1 EMIN=E E=0.5D0*(E+EMAX) IF(IWR.NE.0) THEN WRITE(6,2000) N,K WRITE(6,2001) NR,NZERO,IOTP,NRT WRITE(6,2002) E WRITE(6,2004) EMIN,EMAX ENDIF C IF(EMAX-EMIN.LT.DELL*DABS(EMIN)) THEN IER=5 WRITE(6,1005) RETURN ENDIF C GO TO 2 ENDIF C **** THE CORRECT NUMBER OF NODES HAS BEEN FOUND. PO=P(IOTP) QO=Q(IOTP) C C **** INWARD SOLUTION C CALL DINW(E,EPS,K,IOTP) IF(IER.GT.0) RETURN C **** MATCHING OF THE OUTWARD AND INWARD SOLUTIONS. FACT=PO/P(IOTP) DO 5 I=IOTP,NRT P(I)=P(I)*FACT 5 Q(I)=Q(I)*FACT QI=Q(IOTP) RLAST=R(ILAST) C **** NORMALIZATION CALL SPLINE(R,P,A,B,C,D,0.0D0,0.0D0,ILAST) CALL INTEG2(R,A,B,C,D,0.0D0,RLAST,SUMP,ILAST) CALL SPLINE(R,Q,A,B,C,D,0.0D0,0.0D0,ILAST) CALL INTEG2(R,A,B,C,D,0.0D0,RLAST,SUMQ,ILAST) SUM=SUMP+SUMQ C **** EIGENVALUE CORRECTION. (179) IF(SUM.LT.1.0D-15) SUM=1.0D0 DE=SL*PO*(QI-QO)/SUM EP=E+DE C IF(DE.LT.0.0D0) THEN ICMAX=1 EMAX=E ENDIF C IF(DE.GT.0.0D0) THEN ICMIN=1 EMIN=E ENDIF C IF(ICMIN.EQ.0.AND.EP.LT.EMIN) THEN EMIN=1.1D0*EMIN IF(EP.LT.EMIN) EP=0.5D0*(E+EMIN) ENDIF IF(ICMIN.EQ.1.AND.EP.LT.EMIN) EP=0.5D0*(E+EMIN) IF(ICMAX.EQ.1.AND.EP.GT.EMAX) EP=0.5D0*(E+EMAX) IF(EP.GT.ESUP) EP=0.5D0*(ESUP+E) C IF(IWR.NE.0) THEN WRITE(6,2000) N,K 2000 FORMAT(/2X,'SUBROUTINE DBOUND. N =',I3,' K =',I3) WRITE(6,2001) NR,NZERO,IOTP,NRT 2001 FORMAT(2X,'NR =',I3,' NZERO =',I3,' IOTP = ',I5, 1' NGP =',I5) WRITE(6,2002) EP 2002 FORMAT(2X,'E NEW = ',1P,D22.15) WRITE(6,2003) E,DE 2003 FORMAT(2X,'E OLD = ',1P,D22.15,' DE = ',D11.4) WRITE(6,2004) EMIN,EMAX 2004 FORMAT(2X,'EMIN = ',1P,D12.5,' EMAX = ',D12.5) ENDIF C IF(EP.GT.ESUP.AND.DABS(E-ESUP).LT.DELL*DABS(ESUP)) THEN IER=5 WRITE(6,1005) 1005 FORMAT(1X,'*** ERROR 5 (IN DBOUND): E OUT OF RANGE.'/5X, 1'(ACCUMULATED ROUND-OFF ERRORS?).') RETURN ENDIF EO=E E=EP IF(DMIN1(DABS(DE),DABS(E-EO)).GT.DABS(E*DELL)) GO TO 2 C **** NORMALIZATION FACT=1.0D0/DSQRT(SUM) DO 6 I=1,ILAST P(I)=P(I)*FACT 6 Q(I)=Q(I)*FACT C C **** EXTRACT THE 'RAD' GRID... C DO 500 I=1,NGP RLOC=RAD(I) CALL FINDI(R,RLOC,NRT,J) IF(RLOC-R(J).LT.R(J+1)-RLOC) THEN PIO(I)=P(J) QIO(I)=Q(J) ELSE PIO(I)=P(J+1) QIO(I)=Q(J+1) ENDIF 500 CONTINUE C IF(DABS(P(NRT)).GT.1.0D-5*DABS(P(IOTP))) THEN IER=3 WRITE(6,1003) ENDIF C RLOC=R(ILAST) CALL FINDI(RAD,RLOC,NGP,ILAST) ILAST=ILAST+1 RETURN END C ************************************************************** C SUBROUTINE SOUTW C ************************************************************** SUBROUTINE SOUTW(E,EPS,L,NR,NZERO,IOTP) C C OUTWARD SOLUTION OF THE SCHRODINGER RADIAL EQUATION FOR A C PIECEWISE CUBIC FIELD. POWER SERIES METHOD. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,NPPG=NDIM+1,NPTG=NDIM+NPPG) COMMON/RADWF/RAD(NDIM),PIO(NDIM),QIO(NDIM),NGP,ILAST,IER COMMON/RGRID/R(NPTG),P(NPTG),Q(NPTG),IND(NPTG),NRT COMMON/VGRID/RG(NPPG),RV(NPPG),VA(NPPG),VB(NPPG),VC(NPPG), 1 VD(NPPG),NVT COMMON/POTEN/RV0,RV1,RV2,RV3 COMMON/SINOUT/PI,QI,PF,QF,RA,RB,RLN,NSTEP,NCHS COMMON/NZT/NZMAX NZERO=0 NZMAX=0 AL=L N1=IOTP-1 C P(1)=0.0D0 Q(1)=0.0D0 DO 2 I=1,N1 RA=R(I) RB=R(I+1) IN=IND(I) RV0=VA(IN) RV1=VB(IN) RV2=VC(IN) RV3=VD(IN) PI=P(I) QI=Q(I) CALL SCH(E,AL,EPS) NZERO=NZERO+NCHS IF(NCHS.GT.NZMAX) NZMAX=NCHS IF(NZERO.GT.NR.AND.E.LT.0.0D0) RETURN P(I+1)=PF Q(I+1)=QF IF(I.EQ.1) GO TO 2 C **** RENORMALIZATION. IF(RLN.GT.0.0D0) THEN FACT=DEXP(-RLN) DO 1 K=1,I P(K)=P(K)*FACT Q(K)=Q(K)*FACT 1 CONTINUE ENDIF 2 CONTINUE RETURN END C ************************************************************** C SUBROUTINE SINW C ************************************************************** SUBROUTINE SINW(E,EPS,L,IOTP) C C INWARD SOLUTION OF THE SCHRODINGER RADIAL EQUATION FOR A C PIECEWISE CUBIC FIELD. POWER SERIES METHOD. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,NPPG=NDIM+1,NPTG=NDIM+NPPG, 1 TRINF=5625.0D0) COMMON/RADWF/RAD(NDIM),PIO(NDIM),QIO(NDIM),NGP,ILAST,IER COMMON/RGRID/R(NPTG),P(NPTG),Q(NPTG),IND(NPTG),NRT COMMON/VGRID/RG(NPPG),RV(NPPG),VA(NPPG),VB(NPPG),VC(NPPG), 1 VD(NPPG),NVT COMMON/POTEN/RV0,RV1,RV2,RV3 COMMON/SINOUT/PI,QI,PF,QF,RA,RB,RLN,NSTEP,NCHS AL=L C **** WKB SOLUTION AT THE OUTER GRID POINT. (168) N=NRT 1 N1=IND(N-1) RN=R(N) RVN=VA(N1)+RN*(VB(N1)+RN*(VC(N1)+RN*VD(N1))) RVNP=VB(N1)+RN*(2.0D0*VC(N1)+RN*3.0D0*VD(N1)) CMU=2.0D0*RN*(RVN-E*RN)+AL*(AL+1) IF(CMU.LE.0.0D0) THEN IER=6 WRITE(6,1006) 1006 FORMAT(1X,'*** ERROR 6 (IN SBOUND): RV(NGP)<<0 OR E>0.', 1/5X,'(CHECK THE INPUT POTENTIAL VALUES).') RETURN ENDIF C **** PRACTICAL INFINITY. (170) IF(CMU.LT.TRINF.OR.N.EQ.IOTP+1) THEN CRAT=1.0D0-RN*(RVN+RN*RVNP-2*E*RN)/(CMU*CMU) P(N)=1.0D0 Q(N)=(0.5D0/RN)*CRAT-DSQRT(CMU)/RN ILAST=N ELSE P(N)=0.0D0 Q(N)=0.0D0 N=N-1 GO TO 1 ENDIF C N1=N-IOTP DO 3 J=1,N1 I=N-J I1=I+1 RA=R(I1) RB=R(I) IN=IND(I) RV0=VA(IN) RV1=VB(IN) RV2=VC(IN) RV3=VD(IN) PI=P(I1) QI=Q(I1) CALL SCH(E,AL,EPS) P(I)=PF Q(I)=QF C **** RENORMALIZATION. IF(RLN.GT.0.0D0) THEN FACT=DEXP(-RLN) DO 2 K=I1,N P(K)=P(K)*FACT Q(K)=Q(K)*FACT 2 CONTINUE ENDIF 3 CONTINUE RETURN END C ************************************************************** C SUBROUTINE DOUTW C ************************************************************** SUBROUTINE DOUTW(E,EPS,K,NR,NZERO,IOTP) C C OUTWARD SOLUTION OF THE DIRAC RADIAL EQUATION FOR A PIECE- C WISE CUBIC FIELD. POWER SERIES METHOD. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,NPPG=NDIM+1,NPTG=NDIM+NPPG) COMMON/RADWF/RAD(NDIM),PIO(NDIM),QIO(NDIM),NGP,ILAST,IER COMMON/RGRID/R(NPTG),P(NPTG),Q(NPTG),IND(NPTG),NRT COMMON/VGRID/RG(NPPG),RV(NPPG),VA(NPPG),VB(NPPG),VC(NPPG), 1 VD(NPPG),NVT COMMON/POTEN/RV0,RV1,RV2,RV3 COMMON/DINOUT/PI,QI,PF,QF,RA,RB,RLN,NSTEP,NCHS COMMON/NZT/NZMAX NZERO=0 NZMAX=0 AK=K N1=IOTP-1 C P(1)=0.0D0 Q(1)=0.0D0 DO 2 I=1,N1 RA=R(I) RB=R(I+1) IN=IND(I) RV0=VA(IN) RV1=VB(IN) RV2=VC(IN) RV3=VD(IN) PI=P(I) QI=Q(I) CALL DIR(E,AK,EPS) NZERO=NZERO+NCHS IF(NCHS.GT.NZMAX) NZMAX=NCHS IF(NZERO.GT.NR.AND.E.LT.0.0D0) RETURN P(I+1)=PF Q(I+1)=QF IF(I.EQ.1) GO TO 2 C **** RENORMALIZATION. IF(RLN.GT.0.0D0) THEN FACT=DEXP(-RLN) DO 1 J=1,I P(J)=P(J)*FACT Q(J)=Q(J)*FACT 1 CONTINUE ENDIF 2 CONTINUE RETURN END C ************************************************************** C SUBROUTINE DINW C ************************************************************** SUBROUTINE DINW(E,EPS,K,IOTP) C C INWARD SOLUTION OF THE DIRAC RADIAL EQUATION FOR A PIECE- C WISE CUBIC FIELD. POWER SERIES METHOD. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NDIM=800,NPPG=NDIM+1,NPTG=NDIM+NPPG, 1 TRINF=5625.0D0,SL=137.036D0) COMMON/RADWF/RAD(NDIM),PIO(NDIM),QIO(NDIM),NGP,ILAST,IER COMMON/RGRID/R(NPTG),P(NPTG),Q(NPTG),IND(NPTG),NRT COMMON/VGRID/RG(NPPG),RV(NPPG),VA(NPPG),VB(NPPG),VC(NPPG), 1 VD(NPPG),NVT COMMON/POTEN/RV0,RV1,RV2,RV3 COMMON/DINOUT/PI,QI,PF,QF,RA,RB,RLN,NSTEP,NCHS C **** ORBITAL ANGULAR MOMENTUM QUANTUM NUMBER. (15) IF(K.LT.0) THEN L=-K-1 ELSE L=K ENDIF AK=K AL=L C **** WKB SOLUTION AT THE OUTER GRID POINT. (177) N=NRT FACT=(E+2.0D0*SL*SL)/(SL*SL) 1 N1=IND(N-1) RN=R(N) RVN=VA(N1)+RN*(VB(N1)+RN*(VC(N1)+RN*VD(N1))) RVNP=VB(N1)+RN*(2.0D0*VC(N1)+RN*3.0D0*VD(N1)) CMU=FACT*RN*(RVN-E*RN)+AL*(AL+1) IF(CMU.LE.0.0D0) THEN IER=6 WRITE(6,1006) 1006 FORMAT(1X,'*** ERROR 6 (IN DBOUND): RV(NGP)<<0 OR E>0.', 1/5X,'(CHECK THE INPUT POTENTIAL VALUES).') RETURN ENDIF C **** PRACTICAL INFINITY. (170) IF(CMU.LT.TRINF.OR.N.EQ.IOTP+1) THEN CRAT=(0.5D0-DSQRT(CMU))/RN-0.25D0*FACT*(RVN+RN*RVNP 1 -2.0D0*RN*E)/CMU P(N)=1.0D0 Q(N)=-SL*(CRAT+AK/RN)/(E+2.0D0*SL*SL) ILAST=N ELSE P(N)=0.0D0 Q(N)=0.0D0 N=N-1 GO TO 1 ENDIF C N1=N-IOTP DO 3 J=1,N1 I=N-J I1=I+1 RA=R(I1) RB=R(I) IN=IND(I) RV0=VA(IN) RV1=VB(IN) RV2=VC(IN) RV3=VD(IN) PI=P(I1) QI=Q(I1) CALL DIR(E,AK,EPS) P(I)=PF Q(I)=QF C **** RENORMALIZATION. IF(RLN.GT.0.0D0) THEN FACT=DEXP(-RLN) DO 2 M=I1,N P(M)=P(M)*FACT Q(M)=Q(M)*FACT 2 CONTINUE ENDIF 3 CONTINUE RETURN END C ************************************************************** C SUBROUTINE SCH C ************************************************************** SUBROUTINE SCH(E,AL,EPS) C C THIS SUBROUTINE SOLVES THE RADIAL SCHRODINGER EQUATION FOR C A CENTRAL FIELD V(R) SUCH THAT C R*V(R) = RV0+RV1*R+RV2*R**2+RV3*R**3 C GIVEN THE BOUNDARY CONDITIONS (I.E. THE VALUE OF THE C RADIAL FUNCTION AND ITS DERIVATIVE) AT RA, THE SOLUTION IN C THE INTERVAL BETWEEN RA AND RB IS GENERATED BY USING A C PIECEWISE POWER SERIES EXPANSION FOR A PARTITION OF THE C INTERVAL, SUITABLY CHOSEN TO ALLOW FAST CONVERGENCE OF THE C SERIES. C C INPUT ARGUMENTS: C E ..................... PARTICLE KINETIC ENERGY C AL .................... ORBITAL ANGULAR MOMENTUM QUANTUM C NUMBER C C INPUT (COMMON POTEN): C RV0, RV1, RV2, RV3 .... POTENTIAL PARAMETERS C C INPUT-OUTPUT (COMMON SINOUT): C RA, RB ................ INTERVAL END POINTS (INPUT) C PI, QI ................ VALUES OF THE RADIAL FUNCTION C AND ITS DERIVATIVE AT RA (INPUT) C PF, QF ................ VALUES OF THE RADIAL FUNCTION C AND ITS DERIVATIVE AT RB (OUTPUT) C RLN ................... DLOG OF THE RE-NORMALIZING FACTOR C NSTEP ................. NUMBER OF STEPS C NCHS .................. NUMBER OF ZEROS OF P(R) IN (RA,RB) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/POTEN/RV0,RV1,RV2,RV3 COMMON/SINOUT/PI,QI,PF,QF,RA,RB,RLN,NSTEP,NCHS COMMON/SSAVE/P0,Q0,P1,Q1,CA(60),R0,R1,NSUM NCHS=0 RLN=0.0D0 C H=RB-RA IF(H.LT.0.0D0) THEN DIRECT=-1.0D0 ELSE DIRECT=1.0D0 ENDIF K=-2 NSTEP=0 C R1=RA P1=PI Q1=QI 1 R0=R1 P0=P1 Q0=Q1 2 IOUT=0 R1=R0+H IF(DIRECT*(RB-R1).LT.DIRECT*1.0D-1*H) THEN R1=RB H=RB-R0 IOUT=1 ENDIF CALL SCH0(E,AL,EPS) C K=K+1 IF(NSUM.GT.15) GO TO 3 IF(K.LT.0) GO TO 4 H=H+H K=0 GO TO 4 3 IF(NSUM.LT.60) GO TO 4 H=0.5D0*H K=-4 GO TO 2 4 NSTEP=NSTEP+1 TST=DABS(P1) IF(TST.GT.1.0D2) THEN C **** RENORMALIZATION. RLN=RLN+DLOG(TST) P1=P1/TST Q1=Q1/TST ENDIF IF(P0*P1.LT.0.0D0.AND.R0.GT.0.0D0) NCHS=NCHS+1 IF(IOUT.EQ.0) GO TO 1 C **** OUTPUT. PF=P1 QF=Q1 RETURN END C ************************************************************** C SUBROUTINE SCH0 C ************************************************************** SUBROUTINE SCH0(E,AL,EPS) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C **** OVERFLOW LEVEL. PARAMETER (OVER=1.0D15) COMMON/POTEN/RV0,RV1,RV2,RV3 COMMON/SSAVE/P0,Q0,P1,Q1,CA(60),R0,R1,NSUM C RVE=RV1-E IF(R0.GT.1.0D-10) GO TO 3 C C **** FIRST INTERVAL. (134-143) C S=AL+1 U0=AL*S U1=2*RV0*R1 U2=2*RVE*R1**2 U3=2*RV2*R1**3 U4=2*RV3*R1**4 UT=U0+U1+U2+U3+U4 C CA(1)=1.0D0 CA(2)=U1*CA(1)/((S+1)*S-U0) CA(3)=(U1*CA(2)+U2*CA(1))/((S+2)*(S+1)-U0) CA(4)=(U1*CA(3)+U2*CA(2)+U3*CA(1)) 1 /((S+3)*(S+2)-U0) C P1=CA(1)+CA(2)+CA(3)+CA(4) Q1=S*CA(1)+(S+1)*CA(2)+(S+2)*CA(3)+(S+3)*CA(4) P2P1=S*(S-1)*CA(1)+(S+1)*S*CA(2)+(S+2)*(S+1)*CA(3) 1 +(S+3)*(S+2)*CA(4) C DO 1 I=5,60 K=I-1 CA(I)=(U1*CA(K)+U2*CA(I-2)+U3*CA(I-3)+U4*CA(I-4)) 1 /((S+K)*(S+K-1)-U0) P1=P1+CA(I) DQ1=(S+K)*CA(I) Q1=Q1+DQ1 P2P1=P2P1+(S+K-1)*DQ1 C **** CHECK OVERFLOW LIMIT. TST=DMAX1(DABS(P1),DABS(Q1),DABS(P2P1)) IF(TST.GT.OVER) THEN NSUM=100 RETURN ENDIF T1=DABS(CA(I)) T2=DABS(R1*R1*(P2P1-UT*P1)) TST1=EPS*DMAX1(DABS(P1),DABS(Q1)/I) TST2=EPS*DMAX1(DABS(P1),DABS(Q1)) IF(T1.LT.TST1.AND.T2.LT.TST2) GO TO 2 1 CONTINUE C **** RENORMALIZATION. (143) 2 NSUM=K+1 Q1=Q1/(P1*R1) P1=1.0D0 RETURN C C **** MIDDLE REGION. (134-139) C 3 CONTINUE H=R1-R0 H2=H*H C RHO=H/R0 U0=AL*(AL+1)+2*R0*(RV0+R0*(RVE+R0*(RV2+R0*RV3))) U1=2*(RV0+R0*(2*RVE+R0*(3*RV2+R0*4*RV3)))*H U2=2*(RVE+R0*(3*RV2+R0*6*RV3))*H2 U3=2*(RV2+R0*4*RV3)*H2*H U4=2*RV3*H2*H2 UT=U0+U1+U2+U3+U4 C CA(1)=P0 CA(2)=Q0*H CA(3)=RHO*RHO*U0*CA(1)/2 CA(4)=RHO*(RHO*(U0*CA(2)+U1*CA(1))-4*CA(3))/6 CAK=(U0-2)*CA(3)+U1*CA(2)+U2*CA(1) CA(5)=RHO*(RHO*CAK-12*CA(4))/12 CAK=(U0-6)*CA(4)+U1*CA(3)+U2*CA(2)+U3*CA(1) CA(6)=RHO*(RHO*CAK-24*CA(5))/20 C P1=CA(1)+CA(2)+CA(3)+CA(4)+CA(5)+CA(6) Q1=CA(2)+2*CA(3)+3*CA(4)+4*CA(5)+5*CA(6) P2P1=2*CA(3)+6*CA(4)+12*CA(5)+20*CA(6) C DO 4 I=7,60 K=I-1 CAK=(U0-(K-2)*(K-3))*CA(I-2)+U1*CA(I-3)+U2*CA(I-4) 1 +U3*CA(I-5)+U4*CA(I-6) CA(I)=RHO*(RHO*CAK-2*(K-1)*(K-2)*CA(K))/(K*(K-1)) P1=P1+CA(I) DQ1=K*CA(I) Q1=Q1+DQ1 P2P1=P2P1+K*(K-1)*CA(I) C **** CHECK OVERFLOW LIMIT. TST=DMAX1(DABS(P1),DABS(Q1),DABS(P2P1)) IF(TST.GT.OVER) THEN NSUM=100 RETURN ENDIF T1=DABS(CA(I)) T2=DABS(R1*R1*P2P1-H2*UT*P1) TST1=EPS*DMAX1(DABS(P1),DABS(Q1)/I) TST2=EPS*DMAX1(DABS(P1),DABS(Q1)) IF(T1.LT.TST1.AND.T2.LT.TST2) GO TO 5 4 CONTINUE C 5 NSUM=K+1 Q1=Q1/H RETURN END C ************************************************************** C SUBROUTINE DIR C ************************************************************** SUBROUTINE DIR(E,AK,EPS) C C THIS SUBROUTINE SOLVES THE RADIAL DIRAC EQUATION FOR A C CENTRAL FIELD V(R) SUCH THAT C R*V(R) = RV0+RV1*R+RV2*R**2+RV3*R**3 C GIVEN THE BOUNDARY CONDITIONS (I.E. THE VALUE OF THE C LARGE AND SMALL RADIAL FUNCTIONS) AT RA, THE SOLUTION IN C THE INTERVAL BETWEEN RA AND RB IS GENERATED BY USING A C PIECEWISE POWER SERIES EXPANSION FOR A PARTITION OF THE C INTERVAL, SUITABLY CHOSEN TO ALLOW FAST CONVERGENCE OF THE C SERIES. C C INPUT ARGUMENTS: C E ..................... PARTICLE KINETIC ENERGY C AK .................... RELATIVISTIC ANGULAR MOMENTUM C QUANTUM NUMBER C C INPUT (COMMON POTEN): C RV0, RV1, RV2, RV3 .... POTENTIAL PARAMETERS C C INPUT-OUTPUT (COMMON DINOUT): C RA, RB ................ INTERVAL END POINTS (INPUT) C PI, QI ................ VALUES OF THE LARGE AND SMALL C RADIAL FUNCTIONS AT RA (INPUT) C PF, QF ................ VALUES OF THE LARGE AND SMALL C RADIAL FUNCTIONS AT RB (OUTPUT) C RLN ................... DLOG OF THE RE-NORMALIZING FACTOR C EPS ................... ESTIMATE OF THE GLOBAL ERROR IN C PF AND QF C NSTEP ................. NUMBER OF STEPS C NCHS .................. NUMBER OF ZEROS OF P(R) IN (RA,RB) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/POTEN/RV0,RV1,RV2,RV3 COMMON/DINOUT/PI,QI,PF,QF,RA,RB,RLN,NSTEP,NCHS COMMON/DSAVE/P0,Q0,P1,Q1,CA(60),CB(60),R0,R1,NSUM NCHS=0 RLN=0.0D0 C H=RB-RA IF(H.LT.0.0D0) THEN DIRECT=-1.0D0 ELSE DIRECT=1.0D0 ENDIF K=-2 NSTEP=0 C R1=RA P1=PI Q1=QI 1 R0=R1 P0=P1 Q0=Q1 2 IOUT=0 R1=R0+H IF(DIRECT*(RB-R1).LT.DIRECT*1.0D-1*H) THEN R1=RB H=RB-R0 IOUT=1 ENDIF CALL DIR0(E,AK,EPS) C K=K+1 IF(NSUM.GT.15) GO TO 3 IF(K.LT.0) GO TO 4 H=H+H K=0 GO TO 4 3 IF(NSUM.LT.60) GO TO 4 H=0.5D0*H K=-4 GO TO 2 4 NSTEP=NSTEP+1 TST=DABS(P1) IF(TST.GT.1.0D2) THEN C **** RENORMALIZATION. RLN=RLN+DLOG(TST) P1=P1/TST Q1=Q1/TST ENDIF IF(P0*P1.LT.0.0D0.AND.R0.GT.0.0D0) NCHS=NCHS+1 IF(IOUT.EQ.0) GO TO 1 C **** OUTPUT. PF=P1 QF=Q1 RETURN END C ************************************************************** C SUBROUTINE DIR0 C ************************************************************** SUBROUTINE DIR0(E,AK,EPS) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C **** SPEED OF LIGHT AND OVERFLOW LEVEL. PARAMETER (SL=137.036D0,OVER=1.0D15) COMMON/POTEN/RV0,RV1,RV2,RV3 COMMON/DSAVE/P0,Q0,P1,Q1,CA(60),CB(60),R0,R1,NSUM C ISIG=1 IF(AK.GT.0.0D0) ISIG=-1 H=R1-R0 H2=H*H RVE=RV1-E C IF(R0.GT.1.0D-10) GO TO 7 C C **** FIRST INTERVAL. (147-164) C U0=RV0/SL U1=RVE*R1/SL U2=RV2*R1**2/SL U3=RV3*R1**3/SL UT=U0+U1+U2+U3 UQ=UT-2*SL*R1 UH=U1-2*SL*R1 IF(DABS(U0).LT.1.0D-10) GO TO 2 C C **** U0.NE.0. (155-159) S=DSQRT(AK*AK-U0*U0) DS=S+S CA(1)=1.0D0 CB(1)=(S+AK)/U0 CAI=U1*CA(1) CBI=UH*CB(1) CA(2)=(U0*CAI+(S+1-AK)*CBI)/(DS+1) CB(2)=(-(S+1+AK)*CAI+U0*CBI)/(DS+1) CAI=U1*CA(2)+U2*CA(1) CBI=UH*CB(2)+U2*CB(1) CA(3)=(U0*CAI+(S+2-AK)*CBI)/(2*(DS+2)) CB(3)=(-(S+2+AK)*CAI+U0*CBI)/(2*(DS+2)) P1=CA(1)+CA(2)+CA(3) PP1=S*CA(1)+(S+1)*CA(2)+(S+2)*CA(3) Q1=CB(1)+CB(2)+CB(3) QP1=S*CB(1)+(S+1)*CB(2)+(S+2)*CB(3) C DO 1 I=4,60 K=I-1 CAI=U1*CA(K)+U2*CA(I-2)+U3*CA(I-3) CBI=UH*CB(K)+U2*CB(I-2)+U3*CB(I-3) CA(I)=(U0*CAI+(S+K-AK)*CBI)/(K*(DS+K)) CB(I)=(-(S+K+AK)*CAI+U0*CBI)/(K*(DS+K)) P1=P1+CA(I) PP1=PP1+(S+K)*CA(I) Q1=Q1+CB(I) QP1=QP1+(S+K)*CB(I) C **** CHECK OVERFLOW LIMIT. TST=DMAX1(DABS(P1),DABS(Q1),DABS(PP1),DABS(QP1)) IF(TST.GT.OVER) THEN NSUM=100 RETURN ENDIF T1A=DABS(R1*PP1+H*(AK*P1-UQ*Q1)) T1B=DABS(R1*QP1-H*(AK*Q1-UT*P1)) T1=DMAX1(T1A,T1B) T2=DMAX1(DABS(CA(I)),DABS(CB(I))) TST=EPS*DMAX1(DABS(P1),DABS(Q1)) IF(T1.LT.TST.AND.T2.LT.TST) GO TO 6 1 CONTINUE GO TO 6 C C **** U0.EQ.0 AND SIG=1. (160,161) 2 IF(ISIG.LT.0) GO TO 4 S=DABS(AK) DS1=S+S+1 CA(1)=1.0D0 CB(1)=-U1*CA(1)/DS1 CA(2)=0.0D0 CB(2)=-U2*CA(1)/(DS1+1) CA(3)=UH*CB(1)/2 CB(3)=-(U1*CA(3)+U3*CA(1))/(DS1+2) CA(4)=(UH*CB(2)+U2*CB(1))/3 CB(4)=-(U1*CA(4)+U2*CA(3))/(DS1+3) P1=CA(1)+CA(2)+CA(3)+CA(4) PP1=S*CA(1)+(S+1)*CA(2)+(S+2)*CA(3)+(S+3)*CA(4) Q1=CB(1)+CB(2)+CB(3)+CB(4) QP1=(S+1)*CB(1)+(S+2)*CB(2)+(S+3)*CB(3) C DO 3 I=5,60 K=I-1 CA(I)=(UH*CB(I-2)+U2*CB(I-3)+U3*CB(I-4))/K CB(I)=-(U1*CA(I)+U2*CA(K)+U3*CA(I-2))/(DS1+K) P1=P1+CA(I) PP1=PP1+(S+K)*CA(I) Q1=Q1+CB(I) QP1=QP1+(S+I)*CB(I) C **** CHECK OVERFLOW LIMIT. TST=DMAX1(DABS(P1),DABS(Q1),DABS(PP1),DABS(QP1)) IF(TST.GT.OVER) THEN NSUM=100 RETURN ENDIF T1A=DABS(R1*PP1+H*(AK*P1-UQ*Q1)) T1B=DABS(R1*QP1-H*(AK*Q1-UT*P1)) T1=DMAX1(T1A,T1B) T2=DMAX1(DABS(CA(I)),DABS(CB(I))) TST=EPS*DMAX1(DABS(P1),DABS(Q1)) IF(T1.LT.TST.AND.T2.LT.TST) GO TO 6 3 CONTINUE GO TO 6 C C **** U0.EQ.0 AND SIG=-1. (162,163) 4 S=DABS(AK)+1 DS1=S+DABS(AK) CB(1)=1.0D0 CA(1)=UH*CB(1)/DS1 CB(2)=0.0D0 CA(2)=U2*CB(1)/(DS1+1) CB(3)=-U1*CA(1)/2 CA(3)=(UH*CB(3)+U3*CB(1))/(DS1+2) CB(4)=-(U1*CA(2)+U2*CA(1))/3 CA(4)=(UH*CB(4)+U2*CB(3))/(DS1+3) P1=CA(1)+CA(2)+CA(3)+CA(4) PP1=S*CA(1)+(S+1)*CA(2)+(S+2)*CA(3)+(S+3)*CA(4) Q1=CB(1)+CB(2)+CB(3)+CB(4) QP1=(S-1)*CB(1)+S*CB(2)+(S+1)*CB(3) C DO 5 I=5,60 K=I-1 CB(I)=-(U1*CA(I-2)+U2*CA(I-3)+U3*CA(I-4))/K CA(I)=(UH*CB(I)+U2*CB(K)+U3*CB(I-2))/(DS1+K) P1=P1+CA(I) PP1=PP1+(S+K)*CA(I) Q1=Q1+CB(I) QP1=QP1+(S+K-1)*CB(I) C **** CHECK OVERFLOW LIMIT. TST=DMAX1(DABS(P1),DABS(Q1),DABS(PP1),DABS(QP1)) IF(TST.GT.OVER) THEN NSUM=100 RETURN ENDIF T1A=DABS(R1*PP1+H*(AK*P1-UQ*Q1)) T1B=DABS(R1*QP1-H*(AK*Q1-UT*P1)) T1=DMAX1(T1A,T1B) T2=DMAX1(DABS(CA(I)),DABS(CB(I))) TST=EPS*DMAX1(DABS(P1),DABS(Q1)) IF(T1.LT.TST.AND.T2.LT.TST) GO TO 6 5 CONTINUE C **** RENORMALIZATION. (164) 6 NSUM=K+1 Q1=Q1/P1 P1=1.0D0 RETURN C C **** MIDDLE REGION. (148-152) C 7 CONTINUE RHO=H/R0 U0=(RV0+R0*(RVE+R0*(RV2+R0*RV3)))/SL U1=(RVE+R0*(2*RV2+R0*3*RV3))*H/SL U2=(RV2+R0*3*RV3)*H2/SL U3=RV3*H*H2/SL UB=U0-2*SL*R0 UH=U1-2*SL*H UT=U0+U1+U2+U3 UQ=UT-2*SL*R1 C CA(1)=P0 CB(1)=Q0 CA(2)=RHO*(-AK*CA(1)+UB*CB(1)) CB(2)=-RHO*(U0*CA(1)-AK*CB(1)) CA(3)=RHO*(-(1+AK)*CA(2)+UB*CB(2)+UH*CB(1))/2 CB(3)=-RHO*(U0*CA(2)+(1-AK)*CB(2)+U1*CA(1))/2 CAI=-(2+AK)*CA(3)+UB*CB(3)+UH*CB(2)+U2*CB(1) CBI=U0*CA(3)+(2-AK)*CB(3)+U1*CA(2)+U2*CA(1) CA(4)=RHO*CAI/3 CB(4)=-RHO*CBI/3 C P1=CA(1)+CA(2)+CA(3)+CA(4) PP1=CA(2)+2*CA(3)+3*CA(4) Q1=CB(1)+CB(2)+CB(3)+CB(4) QP1=CB(2)+2*CB(3)+3*CB(4) C DO 9 I=5,60 K=I-1 CAI=-(K-1+AK)*CA(K)+UB*CB(K)+UH*CB(I-2) 1 +U2*CB(I-3)+U3*CB(I-4) CBI=U0*CA(K)+(K-1-AK)*CB(K)+U1*CA(I-2) 1 +U2*CA(I-3)+U3*CA(I-4) CA(I)=RHO*CAI/K CB(I)=-RHO*CBI/K P1=P1+CA(I) PP1=PP1+K*CA(I) Q1=Q1+CB(I) QP1=QP1+K*CB(I) C **** CHECK OVERFLOW LIMIT. TST=DMAX1(DABS(P1),DABS(Q1),DABS(PP1),DABS(QP1)) IF(TST.GT.OVER) THEN NSUM=100 RETURN ENDIF T1A=DABS(R1*PP1+H*(AK*P1-UQ*Q1)) T1B=DABS(R1*QP1-H*(AK*Q1-UT*P1)) T1=DMAX1(T1A,T1B) T2=DMAX1(DABS(CA(I)),DABS(CB(I))) TST=EPS*DMAX1(DABS(P1),DABS(Q1)) IF(T1.LT.TST.AND.T2.LT.TST) GO TO 10 9 CONTINUE C 10 NSUM=K+1 RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCC CUBIC SPLINE INTERPOLATION CCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ************************************************************** C SUBROUTINE SPLINE C ************************************************************** SUBROUTINE SPLINE(X,Y,A,B,C,D,S1,SN,N) C C CUBIC SPLINE INTERPOLATION BETWEEN TABULATED DATA. C INPUT: C X(I) (I=1, ...,N) ........ GRID POINTS. C (THE X VALUES MUST BE IN INCREASING ORDER). C Y(I) (I=1, ...,N) ........ CORRESPONDING FUNCTION VALUES. C S1,SN ..... SECOND DERIVATIVES AT X(1) AND X(N). C (THE NATURAL SPLINE CORRESPONDS TO TAKING S1=SN=0). C N ........................ NUMBER OF GRID POINTS. C THE INTERPOLATING POLYNOMIAL IN THE I-TH INTERVAL, FROM C X(I) TO X(I+1), IS C PI(X) = A(I)+X*(B(I)+X*(C(I)+X*D(I))) C OUTPUT: C A(I),B(I),C(I),D(I) ...... SPLINE COEFFICIENTS. C C REF.: M.J. MARON, 'NUMERICAL ANALYSIS: A PRACTICAL C APPROACH', MACMILLAN PUBL. CO., NEW YORK 1982. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),Y(N),A(N),B(N),C(N),D(N) IF(N.LT.4) THEN WRITE(6,10) N 10 FORMAT(5X,'SPLINE INTERPOLATION CANNOT BE PERFORMED WITH', 1I4,' POINTS. STOP.') STOP ENDIF N1=N-1 N2=N-2 C **** AUXILIARY ARRAYS H(=A) AND DELTA(=D). DO 1 I=1,N1 IF(X(I+1)-X(I).LT.1.0D-10) THEN WRITE(6,11) 11 FORMAT(5X,'SPLINE X VALUES NOT IN INCREASING ORDER. STOP.') STOP ENDIF A(I)=X(I+1)-X(I) 1 D(I)=(Y(I+1)-Y(I))/A(I) C **** SYMMETRIC COEFFICIENT MATRIX (AUGMENTED). DO 2 I=1,N2 B(I)=2.0D0*(A(I)+A(I+1)) K=N1-I+1 2 D(K)=6.0D0*(D(K)-D(K-1)) D(2)=D(2)-A(1)*S1 D(N1)=D(N1)-A(N1)*SN C **** GAUSS SOLUTION OF THE TRIDIAGONAL SYSTEM. DO 3 I=2,N2 R=A(I)/B(I-1) B(I)=B(I)-R*A(I) 3 D(I+1)=D(I+1)-R*D(I) C **** THE SIGMA COEFFICIENTS ARE STORED IN ARRAY D. D(N1)=D(N1)/B(N2) DO 4 I=2,N2 K=N1-I+1 4 D(K)=(D(K)-A(K)*D(K+1))/B(K-1) D(N)=SN C **** SPLINE COEFFICIENTS. SI1=S1 DO 5 I=1,N1 SI=SI1 SI1=D(I+1) H=A(I) HI=1.0D0/H A(I)=(HI/6.0D0)*(SI*X(I+1)**3-SI1*X(I)**3) 1 +HI*(Y(I)*X(I+1)-Y(I+1)*X(I)) 2 +(H/6.0D0)*(SI1*X(I)-SI*X(I+1)) B(I)=(HI/2.0D0)*(SI1*X(I)**2-SI*X(I+1)**2) 1 +HI*(Y(I+1)-Y(I))+(H/6.0D0)*(SI-SI1) C(I)=(HI/2.0D0)*(SI*X(I+1)-SI1*X(I)) 5 D(I)=(HI/6.0D0)*(SI1-SI) RETURN END C ************************************************************** C SUBROUTINE FINDI C ************************************************************** SUBROUTINE FINDI(X,XC,N,I) C C FINDS THE INTERVAL (X(I),X(I+1)) CONTAINING THE VALUE XC. C INPUT: C X(I) (I=1, ...,N) ........ GRID POINTS. C (THE X VALUES MUST BE IN INCREASING ORDER). C XC ....................... POINT TO BE LOCATED. C N ........................ NUMBER OF GRID POINTS. C OUTPUT: C I ........................ INTERVAL INDEX. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N) IF(XC.GT.X(N)) THEN I=N-1 RETURN ENDIF IF(XC.LT.X(1)) THEN I=1 RETURN ENDIF I=1 I1=N 1 IT=(I+I1)/2 IF(XC.GT.X(IT)) I=IT IF(XC.LE.X(IT)) I1=IT IF(I1-I.GT.1) GO TO 1 RETURN END C ************************************************************** C SUBROUTINE INTEG C ************************************************************** SUBROUTINE INTEG(X,A,B,C,D,XL,XU,SUM,N) C C INTEGRAL OF A CUBIC SPLINE FUNCTION. C INPUT: C X(I) (I=1, ...,N) ........ GRID POINTS. C (THE X VALUES MUST BE IN INCREASING ORDER). C A(I),B(I),C(I),D(I) ...... SPLINE COEFFICIENTS. C N ........................ NUMBER OF GRID POINTS. C XL ....................... LOWER LIMIT IN THE INTEGRAL. C XU ....................... UPPER LIMIT IN THE INTEGRAL. C OUTPUT: C SUM ...................... VALUE OF THE INTEGRAL. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),A(N),B(N),C(N),D(N) C **** SET INTEGRATION LIMITS IN INCREASING ORDER. IF(XU.GT.XL) THEN XLL=XL XUU=XU SIGN=1.0D0 ELSE XLL=XU XUU=XL SIGN=-1.0D0 ENDIF C **** CHECK INTEGRAL LIMITS. IWR=0 IF(XLL.LT.X(1).OR.XUU.GT.X(N)) IWR=1 C **** FIND INVOLVED INTERVALS. SUM=0.0D0 CALL FINDI(X,XLL,N,IL) CALL FINDI(X,XUU,N,IU) C **** ONLY A SINGLE INTERVAL INVOLVED. IF(IL.EQ.IU) THEN X1=XLL X2=XUU SUM=X2*(A(IL)+X2*((B(IL)/2)+X2*((C(IL)/3)+X2*D(IL)/4))) 1 -X1*(A(IL)+X1*((B(IL)/2)+X1*((C(IL)/3)+X1*D(IL)/4))) GO TO 2 ENDIF C **** CONTRIBUTIONS FROM DIFFERENT INTERVALS. X1=XLL X2=X(IL+1) SUM=X2*(A(IL)+X2*((B(IL)/2)+X2*((C(IL)/3)+X2*D(IL)/4))) 1 -X1*(A(IL)+X1*((B(IL)/2)+X1*((C(IL)/3)+X1*D(IL)/4))) IL=IL+1 DO 1 I=IL,IU X1=X(I) X2=X(I+1) IF(I.EQ.IU) X2=XUU SUMP=X2*(A(I)+X2*((B(I)/2)+X2*((C(I)/3)+X2*D(I)/4))) 1 -X1*(A(I)+X1*((B(I)/2)+X1*((C(I)/3)+X1*D(I)/4))) 1 SUM=SUM+SUMP 2 SUM=SIGN*SUM C **** INTEGRAL LIMITS OUT OF RANGE. IF(IWR.EQ.1) WRITE(6,10) 10 FORMAT(/'*** WARNING: INTEGRAL LIMITS OUT OF RANGE. ***') RETURN END C ************************************************************** C SUBROUTINE INTEG2 C ************************************************************** SUBROUTINE INTEG2(X,A,B,C,D,XL,XU,SUM,N) C C INTEGRAL OF A SQUARED CUBIC SPLINE FUNCTION. C INPUT: C X(I) (I=1, ...,N) ........ GRID POINTS. C (THE X VALUES MUST BE IN INCREASING ORDER). C A(I),B(I),C(I),D(I) ...... SPLINE COEFFICIENTS. C N ........................ NUMBER OF GRID POINTS. C XL ....................... LOWER LIMIT IN THE INTEGRAL. C XU ....................... UPPER LIMIT IN THE INTEGRAL. C OUTPUT: C SUM ...................... VALUE OF THE INTEGRAL. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),A(N),B(N),C(N),D(N) C **** SET INTEGRATION LIMITS IN INCREASING ORDER. IF(XU.GT.XL) THEN XLL=XL XUU=XU SIGN=1.0D0 ELSE XLL=XU XUU=XL SIGN=-1.0D0 ENDIF C **** CHECK INTEGRAL LIMITS. IWR=0 IF(XLL.LT.X(1).OR.XUU.GT.X(N)) IWR=1 C **** FIND INVOLVED INTERVALS. SUM=0.0D0 CALL FINDI(X,XLL,N,IL) CALL FINDI(X,XUU,N,IU) C **** ONLY A SINGLE INTERVAL INVOLVED. IF(IL.EQ.IU) THEN X1=XLL X2=XUU SUM=X2*(A(IL)*(A(IL)+X2*B(IL)) 1 +X2*X2*(((2*A(IL)*C(IL)+B(IL)**2)/3.0D0) 2 +X2*(((B(IL)*C(IL)+A(IL)*D(IL))/2.0D0) 3 +X2*(((2*B(IL)*D(IL)+C(IL)**2)/5.0D0) 4 +X2*D(IL)*((C(IL)/3.0D0)+X2*D(IL)/7.0D0))))) 5 -X1*(A(IL)*(A(IL)+X1*B(IL)) 6 +X1*X1*(((2*A(IL)*C(IL)+B(IL)**2)/3.0D0) 7 +X1*(((B(IL)*C(IL)+A(IL)*D(IL))/2.0D0) 8 +X1*(((2*B(IL)*D(IL)+C(IL)**2)/5.0D0) 9 +X1*D(IL)*((C(IL)/3.0D0)+X1*D(IL)/7.0D0))))) GO TO 2 ENDIF C **** CONTRIBUTIONS FROM DIFFERENT INTERVALS. X1=XLL X2=X(IL+1) SUM=X2*(A(IL)*(A(IL)+X2*B(IL)) 1 +X2*X2*(((2*A(IL)*C(IL)+B(IL)**2)/3.0D0) 2 +X2*(((B(IL)*C(IL)+A(IL)*D(IL))/2.0D0) 3 +X2*(((2*B(IL)*D(IL)+C(IL)**2)/5.0D0) 4 +X2*D(IL)*((C(IL)/3.0D0)+X2*D(IL)/7.0D0))))) 5 -X1*(A(IL)*(A(IL)+X1*B(IL)) 6 +X1*X1*(((2*A(IL)*C(IL)+B(IL)**2)/3.0D0) 7 +X1*(((B(IL)*C(IL)+A(IL)*D(IL))/2.0D0) 8 +X1*(((2*B(IL)*D(IL)+C(IL)**2)/5.0D0) 9 +X1*D(IL)*((C(IL)/3.0D0)+X1*D(IL)/7.0D0))))) IL=IL+1 DO 1 I=IL,IU X1=X(I) X2=X(I+1) IF(I.EQ.IU) X2=XUU SUMP=X2*(A(I)*(A(I)+X2*B(I)) 1 +X2*X2*(((2*A(I)*C(I)+B(I)**2)/3.0D0) 2 +X2*(((B(I)*C(I)+A(I)*D(I))/2.0D0) 3 +X2*(((2*B(I)*D(I)+C(I)**2)/5.0D0) 4 +X2*D(I)*((C(I)/3.0D0)+X2*D(I)/7.0D0))))) 5 -X1*(A(I)*(A(I)+X1*B(I)) 6 +X1*X1*(((2*A(I)*C(I)+B(I)**2)/3.0D0) 7 +X1*(((B(I)*C(I)+A(I)*D(I))/2.0D0) 8 +X1*(((2*B(I)*D(I)+C(I)**2)/5.0D0) 9 +X1*D(I)*((C(I)/3.0D0)+X1*D(I)/7.0D0))))) IF(SUMP.LT.0.0D0) SUMP=0.0D0 1 SUM=SUM+SUMP 2 SUM=SIGN*SUM C **** INTEGRAL LIMITS OUT OF RANGE. IF(IWR.EQ.1) WRITE(6,10) 10 FORMAT(/'*** WARNING: INTEGRAL LIMITS OUT OF RANGE. ***') RETURN END