C | C Program Vfit-666. By D.W.Knight. Dec 1984. C Solution of Schroedinger Eqn. for Periodic Potential and C Alpha and Omega dependence of Internal Rot. Const. C Reverts to Hamiltonian of Lewis et al., (J Mol Struct (1972), C 12,449) in absence of Distortion Constants. IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/MATRIX/H(100,100),VECT(100,100,2),ESIN(100) COMMON/PARMTR/F0,FN(8),DF(8),VN(12),NF,ND,NV,NP,NBASIS COMMON/FIX/IFIXF0,IFIXF(8),IFIXD(8),IFIXV(12) CHARACTER LABEL(29)*4,TITLE(10)*8,TEXT(50)*8 common/data/text,mmm(100),owr(150),nobs,nang common/info/label common/misc/title c Reading Title (10A8). READ(5,500) TITLE WRITE(6,600) TITLE C Reading Order of Matrix (No. of Basis Funcs.) READ(5,510) NBASIS IF (NBASIS.GT.100) NBASIS=100 IF (NBASIS.LT.12) NBASIS=12 C Reading F0, the Internal Rotation Const. READ(5,520) IFIXF0,F0 NF=0 ND=0 NV=0 DO 10 I=1,8 FN(I)=0.0 IFIXF(I)=0 DF(I)=0.0 10 IFIXD(I)=0 DO 20 I=1,12 VN(I)=0.0 20 IFIXV(I)=0 C Reading Alpha Dependence Parameters Fn. DO 30 I=1,8 READ(5,530) IEND,IFIX,N,P IF(IEND.NE.0)GO TO 40 FN(N)=P IFIXF(N)=IFIX 30 NF=N C Reading Distortion Parameters DFn. 40 DO 50 I=1,8 READ(5,530) IEND,IFIX,N,P IF(IEND.NE.0)GO TO 60 DF(N)=P IFIXD(N)=IFIX 50 ND=N 60 NVDOM=0 VDOM=0.0 C Reading Potential Coeffs. Vn. DO 70 I=1,12 READ(5,530) IEND,IFIX,N,P IF(IEND.NE.0)GO TO 80 VN(N)=P IFIXV(N)=IFIX IF(ABS(P).GT.ABS(VDOM))NVDOM=N IF(NVDOM.EQ.N)VDOM=P 70 NV=N 80 IF(VDOM.EQ.0.0)NVDOM=0 IF(NVDOM.EQ.0)WRITE(6,602) IF(NVDOM.GT.0)WRITE(6,604)NVDOM NP=17+NV C Least squares Fitting Routine. READ(5,510)NCYCLE,ICOV,ICC,NEW IF(NCYCLE.EQ.0)GO TO 90 WRITE(6,610)NBASIS CALL LEASQ(NCYCLE,ICOV,ICC) C Energy Levels, Potential Func, Matrix, Eigenvectors, Derivatives. 90 READ(5,510,END=200)NPRED,IPLOT,IMAT,IVEC,IDER IF(NEW.GT.0)CALL NUFILE(NCYCLE,ICOV,ICC, *NPRED,IPLOT,IMAT,IVEC,IDER) IF(IPLOT.GT.0) CALL VPLOT(IPLOT,NVDOM) IF(NPRED.EQ.0) STOP IF(NPRED.GE.NBASIS)NPRED=NBASIS-1 IF(NCYCLE.GT.0)GO TO 100 WRITE(6,610)NBASIS WRITE(6,620)F0 IF(NF.GT.0)WRITE(6,630)(FN(I),I=1,NF) IF(ND.GT.0)WRITE(6,640)(DF(I),I=1,ND) IF (IPLOT.GT.0) GO TO 100 IF(NV.GT.0)WRITE(6,650)(VN(I),I=1,NV) 100 CALL SETUP(-1.0D0) IF(IMAT.GT.0)CALL MPRINT(0) IV=0 IF(IVEC.GT.0)IV=1 IF(IDER.GT.0)IV=1 CALL HDIAG(NBASIS,IV) DO 110 I=1,NPRED 110 ESIN(I)=H(I,I) CALL SETUP(1.0D0) IF(IMAT.GT.0)CALL MPRINT(1) IF(IVEC.GT.0)IV=2 IF(IDER.GT.0)IV=2 CALL HDIAG(NBASIS,IV) WRITE(6,660) WRITE(6,665)H(1,1) DO 130 I=1,NPRED J=I+1 130 WRITE(6,670)I,ESIN(I),H(J,J) IF(IVEC.GT.0)CALL VPRINT(NBASIS,IVEC) IF(IDER.GT.0)CALL DERPR(NPRED,NVDOM) 200 STOP C 500 FORMAT(10A8) 510 FORMAT(5I3) 520 FORMAT(I2,2X,F16.8) 530 FORMAT(2I1,I2,F16.8) C 600 format(1h ,10a8/1h ,'Torsional Potential Program VFIT') 602 format(1h ,'Free Rotor.') 604 FORMAT(1h ,I3,' Fold Dominated Potential') 610 FORMAT(/1h ,I3,' Basis Functions') 620 FORMAT(/1h ,'F0='F12.8) 630 FORMAT(/1h ,'Coeffs Fn of F(alpha)=F0+SigmaFn*Cos(n*alpha)' */(1h ,5F16.8)) 640 FORMAT(/1h ,'Distortion Coeffs DFn' */(1h ,5F16.8)) 650 FORMAT(/1h ,'Coeffs Vn of V(alpha)=SigmaVn(1-Cos(n*alpha))/2' */(1h ,5F16.8)) 660 FORMAT(/1h ,'Energy Levels'//2h , *'|m| Odd (Sin) Wfn. Even (Cos) Wfn.'/) 665 FORMAT(2h ,'0',16X,F16.8) 670 FORMAT(1h ,I3,2F16.8) END C--------------------------------------------------------------- BLOCK DATA IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/INFO/LABEL(29) CHARACTER*4 LABEL DATA LABEL/' F0',' F1',' F2',' F3',' F4',' F5',' F6', *' F7',' F8',' DF1',' DF2',' DF3',' DF4',' DF5',' DF6',' DF7', *' DF8',' V1',' V2',' V3',' V4',' V5',' V6',' V7',' V8', *' V9',' V10',' V11',' V12'/ END C---------------------------------------------------------------- SUBROUTINE NUFILE(NCYCLE,ICOV,ICC,NPRED,IPLOT,IMAT,IVEC,IDER) C Outputs an Updated Version of the Input File to Device 09. IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/INFO/L0,LF(8),LD(8),LV(12) CHARACTER*4 L0,LF,LD,LV COMMON/MISC/TITLE(10) CHARACTER*8 TITLE COMMON/PARMTR/F0,FN(8),DF(8),VN(12),NF,ND,NV,NP,NBASIS COMMON/FIX/IF0,IFF(8),IFD(8),IFV(12) COMMON/DATA/TEXT(50),MU(50),ML(50),OBS(50),W(50),RES(50), *NOBS,NANG CHARACTER*8 TEXT IF(NCYCLE.EQ.0) RETURN WRITE(9,900)TITLE WRITE(9,910)NBASIS WRITE(9,920)IF0,F0,L0 IF(NF.EQ.0)GO TO 20 DO 10 I=1,NF 10 WRITE(9,925)IFF(I),I,FN(I),LF(I) 20 IF(NF.LT.8)WRITE(9,930) IF(ND.EQ.0)GO TO 40 DO 30 I=1,ND 30 WRITE(9,925)IFD(I),I,DF(I),LD(I) 40 IF(ND.LT.8)WRITE(9,930) IF(NV.EQ.0)GO TO 60 DO 50 I=1,NV 50 WRITE(9,925)IFV(I),I,VN(I),LV(I) 60 IF(NV.LT.12)WRITE(9,930) NEW=1 WRITE(9,910)NCYCLE,ICOV,ICC,NEW IF(NANG.EQ.0)GO TO 68 DO 62 I=1,NANG IF(W(I).EQ.0.0)W(I)=1.0D-6 U=1.0/SQRT(W(I)) 62 WRITE(9,935)TEXT(I),OBS(I),U,RES(I) 68 WRITE(9,950) IE1=NANG+1 DO 70 I=IE1,NOBS U=1.0/SQRT( W(I) ) 70 WRITE(9,940)TEXT(I),MU(I),ML(I),OBS(I),U,RES(I) WRITE(9,950) WRITE(9,910)NPRED,IPLOT,IMAT,IVEC,IDER RETURN C 900 FORMAT(10A8) 910 FORMAT(5I3) 920 FORMAT(I2,2X,F16.8,2X,A4) 925 FORMAT(2I2,F16.8,2X,A4) 930 FORMAT('1') 935 FORMAT(A8,8X,3F16.8) 940 FORMAT(A8,1h ,I3,1h ,I3,3F16.8) 950 FORMAT(1h ) END C---------------------------------------------------------------- SUBROUTINE VPLOT(ISTEP,NDOM) Calculates Points for Plotting the Potential Func. V(alpha). IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/PARMTR/F(17),V(12),NF,ND,NV,NP,NBASIS IF(NDOM.EQ.0)RETURN WRITE(6,600) WRITE(6,610)(I,I=1,NV) WRITE(6,620)(V(I),I=1,NV) WRITE(6,630) IF(ISTEP.GT.360) ISTEP=30 IALFA=0 DO 20 I=1,360 RALFA=FLOAT(IALFA)*1.74532925199D-02 VALFA=0.0 DO 10 N=1,NV 10 VALFA=VALFA + V(N)*(1.0-DCOS(FLOAT(N)*RALFA) )/2.0 WRITE(6,640)IALFA,VALFA IALFA=IALFA+ISTEP IF(IALFA.GE.360) RETURN 20 CONTINUE C 600 FORMAT(/1h ,'Potential Function ', *'V(alpha)=SigmaVn(1-Cos(n*alpha))/2') 610 FORMAT(1h ,'n'5X,I2,11(8X,I2)) 620 FORMAT(1h ,'Vn',12F10.4) 630 FORMAT(/1h ,'alpha',6X,'V(alpha)') 640 FORMAT(1h ,I4,F15.5) END C-------------------------------------------------------------- SUBROUTINE MPRINT(L) C Prints Out Initial Matrix Elements. IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/MATRIX/H(100,100),V(100,100,2),E(100) IF (L.EQ.0) WRITE(6,600) IF (L.EQ.1) WRITE(6,610) DO 20 I=1,12 20 WRITE(6,620)(H(I,J),J=1,12) RETURN 600 FORMAT(/1h ,'Hamiltonian Matrix in the Sin(m*alpha) Basis. ', *'Top Left is (1|H|1).'/) 610 FORMAT(/1h ,'Hamiltonian Matrix in the Cos(m*alpha) Basis. ', *'Top Left is (0|H|0).'/) 620 FORMAT(1h ,12F10.3) END C---------------------------------------------------------------- SUBROUTINE VPRINT(NB,IV) C Prints Eigenvectors in Blocks of 12. IV Specifies Start of Block. IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/MATRIX/H(100,100),V(100,100,2),E(100) DIMENSION IDENT(12) IF(IV.GE.NB) IV=NB-11 IE=IV+11 WRITE(6,600) WRITE(6,630)(I,I=IV,IE) DO 10 I=1,NB 10 WRITE(6,620)(V(I,K,1),K=IV,IE) WRITE(6,610) DO 15 I=IV,IE 15 IDENT(I)=I-1 WRITE(6,630)(IDENT(I),I=IV,IE) DO 20 I=1,NB 20 WRITE(6,620)(V(I,K,2),K=IV,IE) RETURN 600 FORMAT(/1h ,'Eigenvectors Transforming From the Sin Basis'/) 610 FORMAT(/1h ,'Eigenvectors Transforming From the Cos Basis'/) 620 FORMAT(1h ,12F10.6) 630 FORMAT(2X,I3,11(7X,I3)) END C----------------------------------------------------------------- SUBROUTINE DERPR(NPRED,NVDOM) C Prints Out Energy Level Derivatives w.r.t Parameters. IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/PARMTR/P(29),NF,ND,NV,MAXV,NBASIS COMMON/FIX/IFIX(29) COMMON/DERIVS/DER(50,29) CHARACTER*4 DL(29) DATA DL/' dE/dF0 ',' dE/dF1 ',' dE/dF2 ',' dE/dF3 ',' dE/dF4 ', *' dE/dF5 ',' dE/dF6 ',' dE/dF7 ',' dE/dF8 ',' dE/dDF1', *' dE/dDF2',' dE/dDF3',' dE/dDF4',' dE/dDF5',' dE/dDF6', *' dE/dDF7',' dE/dDF8',' dE/dV1 ',' dE/dV2 ',' dE/dV3 ', *' dE/dV4 ',' dE/dV5 ',' dE/dV6 ',' dE/dV7 ',' dE/dV8 ', *' dE/dV9 ',' dE/dV10',' dE/dV11',' dE/dV12'/ IFIX(1)=1 DO 10 I=2,29 10 IF(P(I).NE.0.0)IFIX(I)=1 IF(NV.EQ.0)GO TO 30 DO 20 I=18,MAXV IF(IFIX(I).GT.0)MINV=I IF(IFIX(I).GT.0)GO TO 30 20 CONTINUE 30 IF(ND.EQ.0)GO TO 50 MAXD=ND+9 DO 40 I=10,MAXD IF(IFIX(I).GT.0)MIND=I IF(IFIX(I).GT.0)GO TO 50 40 CONTINUE 50 IF(NF.EQ.0)GO TO 65 MAXF=NF+1 DO 60 I=2,MAXF IF(IFIX(I).GT.0)MINF=I IF(IFIX(I).GT.0)GO TO 65 60 CONTINUE 65 IF(NV.GT.0)GO TO 68 MAXV=18 MINV=18 NV=1 68 IF(ND.GT.0)GO TO 70 MAXD=10 MIND=10 ND=1 70 IF(NF.GT.0)GO TO 71 MAXF=NVDOM+1 IF(NVDOM.EQ.0)MAXF=2 IF(MAXF.GT.9)MAXF=2 MINF=MAXF NF=NVDOM IF(NF.GT.8)NF=1 IF(NVDOM.EQ.0)NF=1 71 DO 72 I=MINF,MAXF 72 IFIX(I)=1 DO 73 I=MIND,MAXD 73 IFIX(I)=1 DO 74 I=MINV,MAXV 74 IFIX(I)=1 WRITE(6,600) ISIGN=-1 ISHIFT=0 DO 90 IW=1,2 IF(IW.EQ.1)WRITE(6,610) IF(IW.EQ.2)WRITE(6,620) WRITE(6,630)DL(1),(DL(J),J=MINF,MAXF),(DL(K),K=MIND,MAXD) *,(DL(L),L=MINV,MAXV) DO 80 I=1,NPRED M=I*ISIGN-ISHIFT DO 75 II=1,MAXV 75 DER(1,II)=0.0 CALL DERCAL(1,M,1.0D0) 80 WRITE(6,640)M,DER(1,1),(DER(1,J),J=MINF,MAXF),(DER(1,K), *K=MIND,MAXD),(DER(1,L),L=MINV,MAXV) ISIGN=1 90 ISHIFT=1 RETURN 600 FORMAT(/1h ,'Energy Level Derivatives w.r.t Parameters.') 610 FORMAT(/1h ,'Sin Block.') 620 FORMAT(/1h ,'Cos Block.') 630 FORMAT(/4X,'m',12(2X,A8)/5X,12(2X,A8)/5X,5(2X,A8)) 640 FORMAT(1h ,I4,12F10.3/5X,12F10.3/5X,5F10.3) END C ----------------------------------------------------------------- SUBROUTINE LEASQ(NCYCLE,ICOV,ICC) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/MATRIX/H(100,100),VECT(100,100,2),ESIN(100) COMMON/PARMTR/P(29),NF,ND,NV,NP,NBASIS COMMON/FIX/IFIX(29) COMMON/DERIVS/DER(50,29) COMMON/NORM/A(29,29) DIMENSION Y(29),DP(29) COMMON/DATA/TEXT(50),MU(50),ML(50),OBS(50),WEIGHT(50),RES(50), *NOBS,NANG COMMON/INFO/LABEL(29) dimension alfa(2) CHARACTER LABEL*4,ALFA*8,TEXT*8 DATA ALFA/'alphamin','alphamax'/ IF(ICC.EQ.0)ICC=10 JCC=ICC CC=FLOAT(JCC*JCC*JCC) C Writing out the Starting Parameters. WRITE(6,630) DO 5 I=1,NP IF(P(I).EQ.0.0)GO TO 5 IF(IFIX(I).EQ.0)WRITE(6,640)LABEL(I),P(I) IF(IFIX(I).GT.0)WRITE(6,650)LABEL(I),P(I) 5 CONTINUE C Reading Input Data, Deciding Which Wang Blocks are Needed, Calculating Weights and Degrees of Freedom. NEFF=0 DO 6 I=1,10 READ(5,510)TEXT(I),OBS(I),U IF(OBS(I).EQ.0.0)GO TO 7 MU(I)=1 WEIGHT(I)=1.0 IF(U.GE.1.0D-6)WEIGHT(I)=1/(U*U) 6 IF(WEIGHT(I).GT.1.0D-4)NEFF=NEFF+1 7 NANG=I-1 IE1=I ICOS=0 ISIN=0 MMAX=0 DO 10 I=IE1,50 READ(5,500,END=20)TEXT(I),MU(I),ML(I),OBS(I),U MMU=ABS(MU(I)) MML=ABS(ML(I)) IF((MMU+MML).EQ.0)GO TO 20 IF(MU(I).GE.0)ICOS=1 IF(ML(I).GE.0)ICOS=1 IF(MU(I).LT.0)ISIN=1 IF(ML(I).LT.0)ISIN=1 IF(MMU.GT.MMAX)MMAX=MMU IF(MML.GT.MMAX)MMAX=MML WEIGHT(I)=1.0 IF(U.GE.1.0D-8)WEIGHT(I)=1/(U*U) 10 IF(WEIGHT(I).GT.1.0D-6)NEFF=NEFF+1 20 NOBS=I-1 IF(MMAX.EQ.0)RETURN IF(MMAX.GE.NBASIS)NBASIS=MMAX+1 IF(NBASIS.GT.100)STOP'Illegal Quantum Number in input data.' NVP=0 DO 25 I=1,NP 25 IF(IFIX(I).GT.0)NVP=NVP+1 NDF=NEFF-NVP IF(NDF.LE.0)STOP'More Parameters than Useful Observations.' WRITE(6,605)NDF FNDF=FLOAT(NDF) C Least Squares Fitting Cycle. DO 900 IROUND=1,NCYCLE WRITE(6,600)IROUND Calculating Energies, Residuals and Derivatives. WRITE(6,610) SWRR=0.0 IF(NANG.EQ.0)GO TO 29 DO 27 I=1,NANG DO 26 J=1,NP 26 DER(I,J)=0.0 CALL NEWTON(I,ITS) WRITE(6,615)TEXT(I),ALFA(MU(I)),OBS(I),RES(I),WEIGHT(I),ITS 27 SWRR=SWRR+WEIGHT(I)*RES(I)*RES(I) 29 IF(ISIN.EQ.0)GO TO 40 CALL SETUP(-1.0D0) CALL HDIAG(NBASIS,1) DO 30 I=1,MMAX 30 ESIN(I)=H(I,I) 40 IF(ICOS.EQ.0)GO TO 50 CALL SETUP(1.0D0) CALL HDIAG(NBASIS,2) 50 DO 60 I=IE1,NOBS MMU=ABS(MU(I)) MML=ABS(ML(I)) IF(MU(I).LT.0)EUP=ESIN(MMU) IF(MU(I).GE.0)EUP=H(MMU+1,MMU+1) IF(ML(I).LT.0)ELO=ESIN(MML) IF(ML(I).GE.0)ELO=H(MML+1,MML+1) RES(I)=OBS(I)-EUP+ELO WRITE(6,620)TEXT(I),MU(I),ML(I),OBS(I),RES(I),WEIGHT(I) SWRR=SWRR+WEIGHT(I)*RES(I)*RES(I) DO 55 J=1,NP 55 DER(I,J)=0.0 CALL DERCAL(I,MU(I),1.0D0) CALL DERCAL(I,ML(I),-1.0D0) 60 CONTINUE DO 75 K=1,NP DO 70 J=1,NP 70 A(K,J)=0.0 DP(K)=0.0 Y(K)=0.0 75 CONTINUE SIGMA=SQRT(SWRR/FNDF) WRITE(6,660)SIGMA IF(IROUND.EQ.1)GO TO 78 IF(SIGMA.GT.SIGOLD)WRITE(6,686) IF(SIGMA.LE.SIGOLD)WRITE(6,684) C Setting Up and Solving the Normal Equations. 78 DO 90 I=1,NOBS DO 80 J=1,NP Y(J)=Y(J)+WEIGHT(I)*RES(I)*DER(I,J) DO 80 K=J,NP A(K,J)=A(K,J)+WEIGHT(I)*DER(I,K)*DER(I,J) A(J,K)=A(K,J) 80 CONTINUE 90 CONTINUE DO 100 J=1,NP 100 IF(IFIX(J).EQ.0) A(J,J)=1.0 CALL GAUJDN(NP) DO 110 J=1,NP DO 110 K=1,NP 110 DP(J)=DP(J) + A(J,K)*Y(K) C Parameter Adjustments. WRITE(6,665) IFLAG=0 DS=0.0 DO 120 J=1,NP DS=DS+DP(J)*Y(J) IF(IFIX(J).GT.0) GO TO 115 IF(P(J).EQ.0.0) GO TO 120 WRITE(6,680) LABEL(J),P(J) GO TO 120 115 P(J)=P(J)+DP(J) ESD=SIGMA*SQRT(A(J,J)) WRITE(6,670) LABEL(J),P(J),ESD,DP(J) CCP=ESD/CC IF(ABS(DP(J)).GT.CCP) IFLAG=1 120 CONTINUE SIGNEW=SQRT( (SWRR-DS)/FNDF ) WRITE(6,682) SIGNEW IF(ICOV.GT.1)CALL RHOCAL(NP,NVP) IF(IFLAG.EQ.0) GO TO 950 SIGOLD=SIGMA 900 CONTINUE IF(ICOV.EQ.1)CALL RHOCAL(NP,NVP) RETURN 950 WRITE(6,690) IF(ICOV.EQ.1)CALL RHOCAL(NP,NVP) RETURN C 500 FORMAT(A8,2I4,2F16.8) 510 FORMAT(A8,8X,2F16.8) C 600 FORMAT(/1h ,'Cycle Number',I3/) 605 FORMAT(/1h ,I3,' Degrees of Freedom.') 610 FORMAT(13X,'mu ml',7X,'Obs',14X,'Obs-Calc',9X,'Weight') 615 FORMAT(1h ,A8,2X,A8,2(1h ,F16.8),1h ,G16.5,4X,I3) 620 FORMAT(1h ,A8,3X,I3,1h ,I3,2(1h ,F16.8),1h ,G16.5) 630 FORMAT(/1h ,'Least Squares Fitting Routine' */1h ,'Initial Parameters') 640 FORMAT(1h ,A4,F16.8,' Const') 650 FORMAT(1h ,A4,F16.8) 660 FORMAT(/1h ,'E.S.D. of an Observation =',F16.8,'/Sqrt(Weight)') 665 FORMAT(/2X,'Estimated Parameters',7X,'E.S.D.',7X,'Shift from Old') 670 FORMAT(1h ,A4,3(1h ,F16.8)) 680 FORMAT(1h ,A4,1h ,F16.8,7X,'Const.') 682 FORMAT(/1h ,'New Expected S.D. of Obs.=',F16.8,'/Sqrt(Weight)') 684 FORMAT(1h ,'Converging.') 686 FORMAT(1h ,'Diverging.') 690 FORMAT(/1h ,'Refinement Terminated.') END C-------------------------------------------------------------- SUBROUTINE RHOCAL(NP,NVP) C Prints the Parameter Correlation Coeffs. IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/NORM/A(29,29) COMMON/FIX/IFIX(29) COMMON/INFO/LABEL(29) CHARACTER*4 LABEL DIMENSION IV(29),RTD(29),RHO(29) IF(NVP.EQ.0)RETURN WRITE(6,600) K=1 DO 20 I=1,NVP DO 10 J=K,NP IF(IFIX(J).GT.0) GO TO 15 10 CONTINUE 15 IV(I)=J 20 K=J+1 WRITE(6,610)(LABEL(IV(I)),I=1,NVP) DO 30 I=1,NVP K=IV(I) 30 RTD(I)=SQRT(A(K,K)) DO 50 I=1,NVP DO 40 J=1,I 40 RHO(J)=A(IV(I),IV(J))/(RTD(I)*RTD(J)) 50 WRITE(6,620)(RHO(J),J=1,I) RETURN 600 FORMAT(/1h ,'Correlation Coefficients.') 610 FORMAT(3X,A4,14(4X,A4)/3X,A4,13(4X,A4)) 620 FORMAT(1h ,15F8.4/1h ,14F8.4) END C-------------------------------------------------------------- SUBROUTINE DERCAL(IOBS,M,SIGN) Calculates Derivatives Using the Hellmann-Feynman Theorem and the C fact that Eigenfunc(k)=Sum Over i of Eigenvect(i,k)*Basisfunc(i). IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/MATRIX/H(100,100),VECT(100,100,2),ESIN(100) COMMON/PARMTR/P(29),NF,ND,NV,NP,NBASIS COMMON/FIX/IFIX(29) COMMON/DERIVS/DER(50,29) IVT=2 IF(M.LT.0)IVT=1 WANG=1.0 IF(M.LT.0)WANG=-1.0 ISHIFT=0 IF(M.GE.0)ISHIFT=1 ISTART=1+ISHIFT K=ABS(M)+ISHIFT C dE/dF0 IF(IFIX(1).EQ.0)GO TO 20 DO 10 I=ISTART,NBASIS FM=FLOAT(I-ISHIFT) VV=VECT(I,K,IVT) 10 DER(IOBS,1)=DER(IOBS,1)+SIGN*VV*VV*FM*FM C dE/dFn 20 IF(NF.EQ.0)GO TO 50 NPE=NF+1 DO 40 IP=2,NPE IF(IFIX(IP).EQ.0)GO TO 40 N=IP-1 DO 30 I=ISTART,NBASIS ML=I-ISHIFT FML=FLOAT(ML) DO 30 J=ISTART,NBASIS MR=J-ISHIFT FMR=FLOAT(MR) DELTA1=0.0 IF(ABS(ML-MR).EQ.N)DELTA1=1.0 DELTA2=0.0 IF((ML+MR).EQ.N)DELTA2=1.0 30 DER(IOBS,IP)=DER(IOBS,IP)+SIGN*VECT(I,K,IVT)*VECT(J,K,IVT)* *FML*FMR*(DELTA1-WANG*DELTA2)/2.0 40 CONTINUE C dE/dDFn 50 IF(ND.EQ.0)GO TO 100 IF(IFIX(10).EQ.0)GO TO 70 DO 60 I=ISTART,NBASIS FM=FLOAT(I-ISHIFT) VV=VECT(I,K,IVT) 60 DER(IOBS,10)=DER(IOBS,10)+SIGN*VV*VV*FM*FM*FM*FM 70 IF(IFIX(11).EQ.0)GO TO 90 DO 80 I=ISTART,NBASIS FM=FLOAT(I-ISHIFT) VV=VECT(I,K,IVT) 80 DER(IOBS,11)=DER(IOBS,11)+SIGN*VV*VV*FM*FM*FM*FM*FM*FM 90 NPE=ND+9 C dE/dVn 100 IF(NV.EQ.0)RETURN DO 200 IP=18,NP IF(IFIX(IP).EQ.0)GO TO 200 N=IP-17 V0=VECT(1,K,2) IF(IVT.NE.2)GO TO 120 DER(IOBS,IP)=DER(IOBS,IP)+SIGN*(0.5*V0*V0-V0*VECT(N+1,K,2)* * 0.707106781186548 D0 ) 120 DO 140 I=ISTART,NBASIS ML=I-ISHIFT DO 140 J=ISTART,NBASIS MR=J-ISHIFT DELTA1=0.0 IF(ML.EQ.MR)DELTA1=1.0 DELTA2=0.0 IF((ML+MR).EQ.N)DELTA2=1.0 DELTA3=0.0 IF(ABS(ML-MR).EQ.N)DELTA3=1.0 140 DER(IOBS,IP)=DER(IOBS,IP)+SIGN*VECT(I,K,IVT)*VECT(J,K,IVT)* *(0.5*DELTA1-0.25*(WANG*DELTA2+DELTA3 ) ) 200 CONTINUE RETURN END C------------------------------------------------------------- SUBROUTINE NEWTON(IOBS,ITS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/PARMTR/P(17),V(12),NF,ND,NV,NP,NBASIS COMMON/FIX/IFIX(17),IFV(12) COMMON/DERIVS/DER(50,29) COMMON/DATA/TEXT(50),MU(50),ML(50),OBS(50),WT(50),RES(50), *NOBS,NANG CHARACTER TEXT*8 DATA CONV/1.74532925199D-02/ C Newton-Raphson iteration to find potential turning point. ALAST=OBS(IOBS)*CONV DO 20 I=1,100 D1V=0.0 D2V=0.0 DO 10 J=1,NV FN=FLOAT(J) D1V=D1V+(FN*V(J)*DSIN(FN*ALAST))/2.0 D2V=D2V+(FN*FN*V(J)*DCOS(FN*ALAST))/2.0 10 CONTINUE DELTA=-D1V/D2V ANEXT=ALAST+DELTA IF(ABS(DELTA).LT.1.0D-10)GO TO 30 ALAST=ANEXT 20 CONTINUE WT(IOBS)=0.0 RETURN 30 ITS=I RES(IOBS)=OBS(IOBS)-ANEXT/CONV IF(D2V.GT.0.0)GO TO 40 WT(IOBS)=0.0 MU(IOBS)=2 RETURN Calculation of the dAlpha/dVn coeffs. 40 DO 50 I=1,NV IF(IFV(I).EQ.0)GO TO 50 FN=FLOAT(I) DER(IOBS,I)=-FN*DSIN(FN*ANEXT)/(2.0*D2V*CONV) 50 CONTINUE RETURN END C------------------------------------------------------------ SUBROUTINE SETUP(WANG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/MATRIX/H(100,100),VECT(100,100,2),ESIN(100) COMMON/PARMTR/F0,FN(8),DF(8),VN(12),NF,ND,NV,NP,NBASIS C Wang=-1 for Sin (odd) Basis, Wang=+1 for Cos (even) Basis. DO 10 I=1,NBASIS DO 10 J=1,NBASIS 10 H(I,J)=0.0 SV=0.0 IF(NV.LT.1)GO TO 25 DO 20 I=1,NV 20 SV=SV+VN(I) SV=SV/2.0 25 ISHIFT=0 IF(WANG.GT.0.0)ISHIFT=1 ISTART=1+ISHIFT C For cos basis H(1,1)=(0|H|0). For sin basis H(1,1)=(1|H|1). DO 30 I=1,NBASIS M=I-ISHIFT FM=FLOAT(M) 30 H(I,I)=F0*FM*FM+SV IF(NV.LT.1)GO TO 80 IF(WANG.LT.0.0)GO TO 50 RT8=SQRT(8.0) MAX=NV+1 DO 40 I=2,MAX M=I-1 40 H(1,I)=-VN(M)/RT8 50 DO 65 I=ISTART,NBASIS ML=I-ISHIFT K=I+1 DO 60 J=K,NBASIS MR=J-ISHIFT IF((MR-ML).GT.NV)GO TO 62 H(I,J)=H(I,J)-VN(MR-ML)/4.0 60 CONTINUE 62 CONTINUE 65 CONTINUE DO 75 I=ISTART,NBASIS ML=I-ISHIFT DO 70 J=I,NBASIS MR=J-ISHIFT IF((ML+MR).GT.NV)GO TO 72 H(I,J)=H(I,J)-WANG*VN(ML+MR)/4.0 70 CONTINUE 72 CONTINUE 75 CONTINUE 80 IF(NF.LT.1)GO TO 110 DO 95 I=ISTART,NBASIS ML=I-ISHIFT FML=FLOAT(ML) K=I+1 DO 90 J=K,NBASIS MR=J-ISHIFT IF((MR-ML).GT.NF)GO TO 92 FMR=FLOAT(MR) H(I,J)=H(I,J)+FML*FMR*FN(MR-ML)/2.0 90 CONTINUE 92 CONTINUE 95 CONTINUE DO 105 I=ISTART,NBASIS ML=I-ISHIFT FML=FLOAT(ML) DO 100 J=I,NBASIS MR=J-ISHIFT IF((ML+MR).GT.NF)GO TO 102 FMR=FLOAT(MR) H(I,J)=H(I,J)-WANG*FML*FMR*FN(ML+MR)/2.0 100 CONTINUE 102 CONTINUE 105 CONTINUE 110 IF(ND.LT.1)RETURN DO 120 I=ISTART,NBASIS FM=FLOAT(I-ISHIFT) FM4=FM*FM*FM*FM FM6=FM4*FM*FM 120 H(I,I)=H(I,I)+DF(1)*FM4 +DF(2)*FM6 RETURN END C----------------------------------------------------------- SUBROUTINE GAUJDN(N) C GAUSS-JORDAN ALGORITHM FOR INVERSION OF A POSITIVE DEFINITE MATRIX. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION H(29) COMMON/NORM/A(29,29) IF (N.LE.1) GO TO 8 K=N DO 5 KK=1,N P=A(1,1) IF (P.LE.0.0)STOP'Normal matrix fails to invert' DO 3 I=2,N Q=A(I,1) QM=-Q IF (I.GT.K) QM=Q H(I)=QM/P DO 2 J=2,I 2 A(I-1,J-1)=A(I,J)+Q*H(J) 3 CONTINUE A(N,N)=1.0/P DO 4 I=2,N 4 A(N,I-1)=H(I) 5 K=K-1 DO 6 J=1,N DO 6 K=1,J 6 A(K,J)=A(J,K) RETURN 8 A(1,1)=1.0/A(1,1) RETURN END C--------------------------------------------------------------- SUBROUTINE HDIAG(N,IVECT) C Jacobi Diagonalisation for Symmetric Matrices. Only the Upper C Triangle is Used. N=Order of Matrix, NDM=Max Order of Matrix. C Eigenvalues are not ordered so that Correlation between Basis C Funcs. and Energies is Preserved. C Dimension Statement in Calling Prog. is H(NDM,NDM), C U(NDM,NDM,2). U=Storage For two sets of Eigenvectors Loaded C according to IVECT=1,2. If IVECT=0 Eigenvectors are not calcd. C NDM and NDG=NDM+1 are set by the Data Statement. IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/MATRIX/H(10000),U(10000,2),ESIN(100) DATA NDM,NDG,ERA,ERB/100,101,0.0001,0.00000001/ NR=0 N1=N-1 IF(IVECT.EQ.0) GO TO 97 C INITIALISE EIGENVECTORS IF WANTED KK=1 KD=1 DO 101 I=1,N K=KK DO 100 J=1,N U(K,IVECT)=0 100 K=K+1 U(KD,IVECT)=1.0 KD=KD+NDG 101 KK=KK+NDM 97 IF(N1.LE.0)RETURN C FIND ABSOLUTELY LARGRST ELEMENT OF H HTOP=0. KD=1 DO 103 I=1,N K=KD DO 104 J=I,N Q=DABS(H(K)) IF(HTOP-Q) 105,104,104 105 HTOP=Q 104 K=K+1 103 KD=KD+NDG IF(HTOP.LE.0.0)RETURN CALCULATE THE PIVOT THRESHOLD -- THRSH AVGF=FLOAT(N*N1)*.55 D=0 KS=NDG DO 107 JJ=1,N1 K=KS DO 108 J=JJ,N1 Q=SNGL(H(K))/HTOP D=D+Q*Q 108 K=K+NDG 107 KS=KS+NDM DSTOP=ERA*D THRSH=SQRT(D/AVGF)*HTOP 110 IFLAG=0 C START A SWEEP KDD=NDM JJT=0 DO 111 JJ=1,N1 K=KDD+1 MS=1 DO 112 J=JJ,N1 HIJ=H(K) AHIJ=DABS(HIJ) CHECK TO SEE IF PIVOT IS ABOVE THRESHOLD IF(AHIJ-THRSH) 113,113,114 114 KI=K-KDD HII=H(KI) KJ=JJ+K HJJ=H(KJ) S=HJJ-HII AS=DABS(S) CHECK TO SEE IF ROTATION IS SIGNIFICANT IF(AHIJ-ERB*AS) 113,113,115 115 IFLAG=1 CHECK FOR ROTATION CLOSE TO 45 DEGREES IF(1.E-10*AHIJ-AS) 116,117,117 117 S=0.707106781186548D+00 C=S IF(HIJ.LT.0) S=-S GO TO 118 CALCULATE ROTATION WHICH IS NOT CLOSE TO 45 DEGREES -- COS=C,SIN=S 116 T=HIJ/S TT=0.25/DSQRT(0.25+T*T) C=DSQRT(0.5+TT) TT=T*TT/C IF(S) 140,141,141 140 S=-C C=TT+TT GO TO 118 141 S=TT+TT CALCULATE NEW ELEMENTS OF H 118 KK=MS+KDD DO 119 M=MS,KI T=H(M) TT=H(KK) H(M)=C*T-S*TT H(KK)=S*T+C*TT 119 KK=KK+1 MI=KI IF(JJT) 122,122,120 120 MSS=KK+1 DO 121 M=MSS,KJ MI=MI+NDM T=H(MI) TT=H(KK) H(MI)=C*T-S*TT H(KK)=S*T+C*TT 121 KK=M 122 H(KJ)=S*HIJ+C*HJJ H(KI)=C*H(KI)-S*(C*HIJ-S*HJJ) M=KJ DO 123 I=J,N1 MI=MI+NDM T=H(MI) TT=H(M) H(MI)=C*T-S*TT H(M)=S*T+C*TT 123 M=M+NDM NR=NR+1 IF(IVECT.EQ.0) GO TO 98 CALCULATE NEW EIGENVECTORS IF NEEDED MSS=MS+N1 DO 125 I=MS,MSS M=I+KDD T=U(I,IVECT) TT=U(M,IVECT) U(I,IVECT)=C*T-S*TT 125 U(M,IVECT)=S*T+C*TT 98 Q=AHIJ/HTOP D=D-Q*Q CALCULATE THRESHOLD FROM SCRATCH IF ROUND-OFF GETTING LARGE IF(D-DSTOP) 126,129,129 126 D=0. MSS=NDG DO 127 KK=1,N1 M=MSS DO 128 I=KK,N1 Q=SNGL(H(M))/HTOP D=D+Q*Q 128 M=M+NDG 127 MSS=MSS+NDM DSTOP=ERA*D CALCULATE NEW THRESHOLD 129 THRSH=SQRT(D/AVGF)*HTOP 113 K=K+NDG 112 MS=MS+NDM JJT=JJ 111 KDD=KDD+NDM C STOP IF THERE WERE NO SIGNIFICANT ROTATIONS AROUND PIVOTS ABOVE C THRESHOLD IF(IFLAG) 110,199,110 199 RETURN end C---------------------------------------------------------------------|