C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C A N H A R M - Energies, eigenvectors and vibrational transitions C for a reduced anharmonic potential C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Original version written by Johan Mjoberg at University College London C ca 1978. C C The listing used as the basis for the present modification was rescued in C 1989 from UCL before the demise of the EUCLID system. C C The modifications of this version relative to the original are fairly C extensive but, with one exception, they mainly concern readability of the C source code, input and output. C The exception is a rewritten relative intensity calculation for realistic C relative intensities. C C-------------------------------------------------------------------------------- C C ANHARM was used to poduce the tables of eigenvalues and eigenvectors C for the reduced quartic-quadratic potential in C APPENDIX II of Z.Kisiel, PhD Thesis, University of London (1980) C C References: C C J.Laane, Applied Spectroscopy 24,73-80 (1970) - tabulations of energies C for quartic-quadratic reduced anharmonic potential and C useful relations C P.J.Mojoberg and J.Almlof, Chem.Phys. 29,201(1978) - examples of use C C C Version 8a.IV.2004 ----- Zbigniew KISIEL ----- C C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C C Input is taken from file ANHARM.INP C Output is written to file ANHARM.OUT C C The data file should be as below, note that this example has been shifted C by three columns for clarity - so that T in THIETANE is in fact in column 1. C For numerical constants the first 26 charcters are ignored on input, the C first two (integer) constants should fit within six characters, the remaining C floating point constants should fit within 15 characters. C C column1 column27 C | | C Line 1: THIETANE: Wieser,Duckett,Kydd, J.Mol.Spectrosc. 51,115(1974) C Line 2: ! C Line 3: Number of basis functions: 50 C Line 4: Number of energy levels: 10 C Line 5: Energy multiplier: 23.33 C Line 6: V2: -6.827 C Line 7: V3: 0.0 C Line 8: V4: 1.0 C Line 9: V6: 0.0 C Line10: Temperature/K: 293.15 C C================================================================================ C PROGRAM ANHARM IMPLICIT REAL*8 (A-H,O-Z) parameter (MAXDIM=100,MAXTR=4*MAXDIM) COMMON E(MAXDIM),EV(MAXDIM,MAXDIM),V2,V3,V4,V6,HNYHF,LN DIMENSION A(4),ITEXT(20) DIMENSION II(maxtr),LL(maxtr),PP(maxtr),QQ(maxtr),RR(maxtr), * DV(maxtr),ipt(maxtr) character cdum COMMON /FRBLK/dv COMMON /LNUMS/IPT C 500 FORMAT(20A4) 600 FORMAT(1x/1x,78(1H-)/1x,20A4/1x,78(1H-) *///' V=',F8.3,'*(',F8.3,' Z^2 +',F8.3,' Z^3 +', * F8.3,' Z^4 +',F8.3,' Z^6 ) CM**(-1)'// * ' NUMBER OF BASIS VECTORS=',I4/13x,'TEMPERATURE=',F8.3,' K'//) 2001 FORMAT(/8x,'ENERGY LEVELS (cm-1):'/) 2000 FORMAT(8x,'E(',I3,')=',F12.3,F13.3) 2500 FORMAT(//8X,'MATRIX ELEMENTS:'//22x, * 'K=1',13X,'K=2',13X,'K=3',13X,'K=4',/) 3000 FORMAT(1X,'(',I2,'/ Z^K /',I2,')',4(4X,1PD12.5)) 3500 FORMAT(//3x,'Delta_v=1'/ * 3x,'TRANSITIONS FREQUENCY Irel1 (%)'/) 4000 FORMAT(7x,I2,'-',I2,5X,F10.3,4X,F8.2) 4500 FORMAT(//3x,'Delta_v>1'/ * 3X,'TRANSITIONS FREQUENCY Irel1 (%) Irel2 (%)'/) 5000 FORMAT(7x,I2,'-',I2,5X,F10.3,4X,F8.2,5X,F8.2) C C open(6,file='anharm.out',status='unknown') write(*,'(1x//)') WRITE(*,3344) WRITE(6,3344) 3344 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | A N H A R M - ENERGIES, EIGENVECTORS AND VIBRATIONAL ', * 'TRANSITIONS',T79,'|'/ * ' | FOR A REDUCED ANHARMONIC POTENTIAL', * T79,'|'/ * ' |',76(1H_),'|'/' version 8a.III.2004',T64,'Zbigniew KISIEL'/ * T64,'(Johan MJOBERG)'/) WRITE(*,601) 601 FORMAT(1x,/' Input from ANHARM.INP'/ * ' Output to ANHARM.OUT'/) C C C...set up input and output and read the data C open(5,file='anharm.inp',status='old',err=501) C C ITEXT - descriptive comment C N - number of basis vectors C M - number of vibrational states to be calculated C HNYHF - multiplier for the reduced potential C V2,V3,V4,V6 - coefficients in the reduced potential C T - temperature (Kelvin) C 1 READ(5,500,END=222) ITEXT read(5,'(a)')cdum READ(5,'(26x,i6)',err=222)N READ(5,'(26x,i6)',err=222)M READ(5,'(26x,F15.0)',err=222)HNYHF READ(5,'(26x,F15.0)',err=222)V2 READ(5,'(26x,F15.0)',err=222)V3 READ(5,'(26x,F15.0)',err=222)V4 READ(5,'(26x,F15.0)',err=222)V6 READ(5,'(26x,F15.0)',err=222)T C if(n.eq.0.or.hnyhf.eq.0.0.or.m.eq.0)goto 222 write(*,'(1x/'' Successful input for:''//1x,20A4/)')ITEXT LN=N WRITE(6,600) ITEXT,HNYHF,V2,V3,V4,V6,N,T C C C...Do the calculation C CALL ANVIB(N) C C C...write the energies C write(6,2001) write(*,2001) DO 50 I=1,M I34=I-1 WRITE(6,2000) I34, E(I) , E(I)-E(1) WRITE(*,2000) I34, E(I) , E(I)-E(1) 50 CONTINUE C C...write the matrix elements C WRITE(6,2500) DO 200 I=1,M DO 100 K=1,4 I35=I-1 A(K)=ELM(K,I35,I35) 100 CONTINUE WRITE(6,3000) I35,I35,(A(J),J=1,4) 200 CONTINUE C C...calculate and write the deltav=1 transitions and intensities C RK=1.0D0/1.4388D0 MM1=M-1 WRITE(6,3500) write(*,3500) c ppimax=0.d0 DO 301 L=1,MM1 I=L-1 ppi=(DEXP(-E(L)/(RK*T))-DEXP(-E(L+1)/(RK*T))) * *(E(L+1)-E(L))*ELM(1,I,L)*ELM(1,I,L) if(ppi.gt.ppimax)ppimax=ppi 301 continue DO 300 L=1,MM1 I=L-1 P=E(L+1)-E(L) c ppi=(DEXP(-E(L)/(RK*T))-DEXP(-E(L+1)/(RK*T)))*P * *ELM(1,I,L)*ELM(1,I,L) c WRITE(6,4000) I,L,P,100.*ppi/ppimax WRITE(*,4000) I,L,P,100.*ppi/ppimax 300 CONTINUE C C...calculate and write the deltav>1 transitions C WRITE(6,4500) MM2=M-2 ntr=0 pp2max=0.0 DO 400 K=1,MM2 I=K-1 KP1=K+1 DO 450 L=KP1,MM1 P=E(L+1)-E(K) Q=(DEXP(-E(K)/(RK*T))-DEXP(-E(L+1)/(RK*T)))*P * *ELM(1,I,L)*ELM(1,I,L) Q=Q/ppimax R=(DEXP(-E(K)/(RK*T))-DEXP(-E(L+1)/(RK*T)))*P * *ELM(2,I,L)*ELM(2,I,L) if(r.gt.pp2max)pp2max=r c c..code inserted to allow sorting of transitions without modifying any of c the original coding ntr=ntr+1 II(ntr)=i LL(ntr)=l PP(ntr)=p QQ(ntr)=q RR(ntr)=r DV(ntr)=real(l-i)+ntr/(2.d0*maxtr) IPT(ntr)=ntr 450 CONTINUE 400 CONTINUE C C...sort according to deltaV and print C call sorth(1,ntr) do 505 k=1,ntr l=ipt(k) if(k.gt.1)then lprev=ipt(k-1) if( (ll(l)-ii(l)).gt.(ll(lprev)-ii(lprev)) )then write(6,'(2x)') endif endif write(6,5000) ii(l),ll(l),PP(l),100.*QQ(l),100.*RR(l)/pp2max 505 continue C write(6,550) 550 format(/ * 3x,'Irel1 normalized to strongest in the Deltav=1 stack'/ * 3x,'Irel2 normalized to strongest in the Deltav=2 stack'// * 79(1H-)/) GO TO 1 c 222 CONTINUE goto 502 C 501 write(*,'(1x//'' ERROR: No ANHARM.INP file''//)') C 502 STOP END C C------------------------------------------------------------------------------ C SUBROUTINE ANVIB(N) C IMPLICIT REAL*8 (A-H,O-Z) parameter (MAXDIM=100,MDIMSQ=MAXDIM*MAXDIM, * MDLIN=(MAXDIM/2)*(MAXDIM+1)) COMMON E(MAXDIM),EV(MAXDIM,MAXDIM),V2,V3,V4,V6,HNYHF,LN DIMENSION H(MAXDIM,MAXDIM) DIMENSION A(MDLIN),R(MDIMSQ) C DO 100 I=1,N DO 100 J=I,N IM1=I-1 JM1=J-1 H(I,J)=HME(0,IM1,JM1)+V2*HME(2,IM1,JM1)+V3*HME(3,IM1,JM1)+ * V4*HME(4,IM1,JM1)+V6*HME(6,IM1,JM1) IF(I.EQ.J) GO TO 100 H(J,I)=H(I,J) 100 CONTINUE C K=0 DO 125 I=1,N DO 125 J=1,N IF (J.GT.I) GO TO 125 K=K+1 A(K)=H(J,I) 125 CONTINUE C CALL EIGEN(A,R,N,0) C II=0 DO 150 I=1,N II=II+I E(I)=A(II) 150 CONTINUE C K=0 DO 175 I=1,N DO 175 J=1,N K=K+1 EV(J,I)=R(K) 175 CONTINUE C DO 200 I=1,LN E(I)=HNYHF*E(I) 200 CONTINUE C RETURN END C C------------------------------------------------------------------------------ C DOUBLE PRECISION FUNCTION ELM(K,NN,MM) c c To evaluate < NN | mu_K | MM > for mu_K=1 c IMPLICIT REAL*8 (A-H,O-Z) parameter (MAXDIM=100) COMMON E(MAXDIM),EV(MAXDIM,MAXDIM),V2,V3,V4,V6,HNYHF,LN C N=NN+1 M=MM+1 ELM=0.0D0 DO 100 I=1,LN DO 100 J=1,LN IM1=I-1 JM1=J-1 ELM=ELM+EV(I,N)*HME(K,IM1,JM1)*EV(J,M) 100 CONTINUE C RETURN END C C------------------------------------------------------------------------------ C DOUBLE PRECISION FUNCTION HME(N,K,L) IMPLICIT REAL*8 (A-H,O-Z) C C Matrix elements C HME=0.0D0 FK=DFLOAT(K) c IF(N.NE.0) GO TO 20 IF(L.NE.K+2) GO TO 11 HME=-0.5D0*DSQRT((FK+1.0D0)*(FK+2.0D0)) GO TO 13 11 IF(L.NE.K) GO TO 12 HME=FK+0.5D0 GO TO 13 12 IF(L.NE.K-2) GO TO 13 HME=-0.5D0*DSQRT(FK*(FK-1.0D0)) 13 RETURN c 20 IF(N.NE.1) GO TO 30 IF(L.NE.K+1) GO TO 21 HME=DSQRT(0.5D0*(FK+1.0D0)) GO TO 22 21 IF(L.NE.K-1) GO TO 22 HME=DSQRT (0.5D0*FK) 22 RETURN c 30 IF(N.NE.2) GO TO 40 IF(L.NE.K+2) GO TO 31 HME=0.5D0*DSQRT((FK+1.0D0)*(FK+2.0D0)) GO TO 33 31 IF(L.NE.K) GO TO 32 HME=FK+0.5D0 GO TO 33 32 IF(L.NE.K-2) GO TO 33 HME=0.5D0*DSQRT(FK*(FK-1.0D0)) 33 RETURN c 40 IF(N.NE.3) GO TO 50 IF(L.NE.K+3) GO TO 41 HME=DSQRT((FK+1.0D0)*(FK+2.0D0)*(FK+3.0D0)/8.0D0) GO TO 44 41 IF(L.NE.K+1) GO TO 42 HME=3.0D0*DSQRT((FK+1.0D0)**3/8.0D0) GO TO 44 42 IF(L.NE.K-1) GO TO 43 HME=3.0D0*DSQRT(FK**3/8.0D0) GO TO 44 43 IF(L.NE.K-3) GO TO 44 HME=DSQRT(FK*(FK-1.0D0)*(FK-2.0D0)/8.0D0) 44 RETURN c 50 IF(N.NE.4) GO TO 70 IF(L.NE.K+4) GO TO 51 HME=0.25D0*DSQRT((FK+1.0D0)*(FK+2.0D0)*(FK+3.0D0)*(FK+4.0D0)) GO TO 55 51 IF(L.NE.K+2) GO TO 52 HME=0.5D0*(2.0D0*FK+3.0D0)*DSQRT((FK+1.0D0)*(FK+2.0D0)) GO TO 55 52 IF(L.NE.K) GO TO 53 HME=0.25D0*(6.0D0*FK**2+6.0D0*FK+3.0D0) GO TO 55 53 IF(L.NE.K-2) GO TO 54 HME=0.5D0*(2.0D0*FK-1.0D0)*DSQRT(FK*(FK-1.0D0)) GO TO 55 54 IF(L.NE.K-4) GO TO 55 HME=0.25D0*DSQRT(FK*(FK-1.0D0)*(FK-2.0D0)*(FK-3.0D0)) 55 RETURN c 70 IF(N.NE.6) RETURN DELTA=1.0D-10 IF(L.NE.K+6) GO TO 71 HME=0.125D0*DSQRT((FK+1.0D0)*(FK+2.0D0)*(FK+3.0D0)*(FK+4.0D0)*(FK+ *5.0D0)*(FK+6.0D0)) GO TO 77 71 IF(L.NE.K+4) GO TO 72 HME=0.375D0*(FK+1.0D0)*(FK+4.0D0)*DSQRT((FK+1.0D0)*(FK+2.0D0)*(FK+ *3.0D0)*(FK+4.0D0)) GO TO 77 72 IF (L.NE.K+2) GO TO 73 HME=0.375D0*DSQRT((FK+1.0D0)*(FK+2.0D0)*(FK+3.0D0)**4)+3.0D0*0.375 *D0*DSQRT(((FK+1.0D0)**3)*((FK+2.0D0)**3))+0.375D0*DSQRT((FK**4)*(F *K+1.0D0)*(FK+2.0D0)) GO TO 77 73 IF(L.NE.K) GO TO 74 HME=0.125D0*(FK+1.0D0)*(FK+2.0D0)*(FK+3.0D0)+1.125D0*(FK+1.0D0)**3 *+1.125D0*FK**3+0.125D0*FK*(FK-1.0D0)*(FK-2.0D0) GO TO 77 74 IF(L.NE.K-2) GO TO 75 HME=0.375D0*DSQRT(FK*(FK-1.0D0)*(FK+1.0D0)**4+DELTA)+1.125D0*DSQRT *((FK**3)*((FK-1.0D0)**3)+DELTA)+0.375D0*DSQRT(FK*(FK-1.0D0)*(FK-2. *0D0)**4+DELTA) GO TO 77 75 IF(L.NE.K-4) GO TO 76 HME=0.375D0*DSQRT((FK**3)*(FK-1.0D0)*(FK-2.0D0)*(FK-3.0D0)+DELTA)+ *0.375D0*DSQRT(FK*(FK-1.0D0)*(FK-2.0D0)*(FK-3.0D0)**3+DELTA) GO TO 77 76 IF(L.NE.K-6) GO TO 77 HME=0.125D0*DSQRT(FK*(FK-1.0D0)*(FK-2.0D0)*(FK-3.0D0)*(FK-4.0D0)*( *FK-5.0D0)+DELTA) C 77 RETURN END C C------------------------------------------------------------------------------ C SUBROUTINE EIGEN(A,R,N,MV) KTH-EIGE C C .................................................................. C C SUBROUTINE EIGEN, NOW MODIFIED FOR DOUBLE PRECISION C C PURPOSE C COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX C C USAGE C CALL EIGEN(A,R,N,MV) C C DESCRIPTION OF PARAMETERS C A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION. C RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF MATRIX A C IN DESCENDING ORDER. C R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, IN C SAME SEQUENCE AS EIGENVALUES). C N - ORDER OF MATRICES A AND R. C MV - INPUT CODE (0=COMPUTE EIGENVALUES AND EIGENVECTORS, C 1=COMPUTE EIGENVALUES ONLY). C C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE C IN C COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION STATEMENT C WHICH FOLLOWS: C IMPLICIT REAL*8 (A-H,O-Z) C C THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO CONTAIN C DOUBLE PRECISION FORTRAN FUNCTIONS. SQRT IN STATEMENTS 40, 68, 75, C AND 78 MUST BE CHANGED TO DSQRT. ABS IN STATEMENT 62 MUST BE C CHANGED TO DABS. THE CONSTANT IN STATEMENT 5 SHOULD BE CHANGED C TO 1.0D-12. C C .................................................................. C C GENERATE IDENTITY MATRIX C parameter (MAXDIM=100,MDIMSQ=MAXDIM*MAXDIM, * MDLIN=(MAXDIM/2)*(MAXDIM+1)) DIMENSION A(MDLIN),R(MDIMSQ) C 5 RANGE = 1.0D-12 NN = N*N NNN = (NN+N)/2 IF (MV.EQ.1) GO TO 25 DO 20 J=1,NN 20 R(J) = 0.0D0 IJ=-N DO 15 J=1,N IJ=IJ+N+1 15 R(IJ) = 1.0D0 C C COMPUTE INITIAL AND FINAL NORMS (ANROM AND ANRMX) C 25 ANORM = 0.0D0 L = 1 DO 36 I=2,N K = I L = L+I-1 LL = L DO 35 J=L,NNN IF (J.LT.LL) GO TO 35 ANORM = ANORM+A(J)*A(J) K = K+1 LL = LL+K 35 CONTINUE 36 CONTINUE IF (ANORM) 165,165,40 40 ANORM = 1.414D0*DSQRT(ANORM) ANRMX = ANORM*RANGE/DFLOAT(N) C C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR C IND = 0 THR = ANORM 45 THR = THR/DFLOAT(N) 50 L = 1 LQ = 0 ILQ = -N 55 M = L+1 ILQ = ILQ+N LQ = LQ+L-1 MQ = LQ LL = L+LQ 60 MQ = MQ+M-1 LM = L+MQ 62 IF (DABS(A(LM))-THR) 130,65,65 65 IND = 1 MM = M+MQ X = 0.5D0*(A(LL)-A(MM)) 68 Y = -A(LM)/DSQRT(A(LM)*A(LM)+X*X) IF (X) 70,75,75 70 Y = -Y 75 SINX = Y/DSQRT(2.0D0*(1.0D0+(DSQRT(DABS(1.0D0-Y*Y))))) SINX2 = SINX*SINX 78 COSX = DSQRT(1.0D0-SINX2) COSX2 = COSX*COSX SINCS = SINX*COSX C C ROTATE L AND M COLUMNS C IMQ = N*(M-1) IQ = 0 DO 125 I=1,N IQ = IQ+I-1 IF (I-L) 85,120,80 80 IF (I-M) 85,120,90 85 IM = I+MQ GO TO 95 90 IM = M+IQ 95 IF (I-L) 100,105,105 100 IL = I+LQ GO TO 110 105 IL = L+IQ 110 X = A(IL)*COSX-A(IM)*SINX A(IM) = A(IL)*SINX+A(IM)*COSX A(IL) = X 115 IF (MV-1) 120,125,120 120 ILR = ILQ+I IMR = IMQ+I X = R(ILR)*COSX-R(IMR)*SINX R(IMR) = R(ILR)*SINX+R(IMR)*COSX R(ILR) = X 125 CONTINUE X = 2.0D0*A(LM)*SINCS Y = A(LL)*COSX2+A(MM)*SINX2-X X = A(LL)*SINX2+A(MM)*COSX2+X A(LM) = (A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL) = Y A(MM) = X C C TESTS FOR COMPLETION C C TEST FOR M = LAST COLUMN C 130 IF (M-N) 135,140,135 135 M = M+1 GO TO 60 C C TEST FOR L = SECOND FROM LAST COLUMN C 140 IF (L-(N-1)) 145,150,145 145 L = L+1 GO TO 55 150 IF (IND-1) 160,155,160 155 IND = 0 GO TO 50 C C COMPARE THRESHOLD WITH FINAL NORM C 160 IF (THR-ANRMX) 165,165,45 C C SORT EIGENVALUES AND EIGENVECTORS C 165 IQ = -N LL = 0 JJQ = -(N+N) DO 185 I=1,N IQ = IQ+N LL = LL+I MM = LL-I JJQ = JJQ+N JQ = JJQ DO 184 J=I,N JQ = JQ+N MM = MM+J IF (A(MM)-A(LL)) 170,184,184 170 X = A(LL) A(LL) = A(MM) A(MM) = X IF (MV-1) 175,184,175 175 DO 180 K=1,N ILR = IQ+K IMR = JQ+K X = R(ILR) R(ILR) = R(IMR) 180 R(IMR) = X 184 CONTINUE 185 CONTINUE C RETURN END C C------------------------------------------------------------------------------ C SUBROUTINE SORTH(NSTART,N) c c This routine is based on the SORT2 'heapsort' routine from Numerical c Recipes and sorts the quantities in vector WK from WK(NSTART) to WK(N) C in ascending order of magnitude and also accordingly rearranges vector C IPT of pointers to original positions of sorted quantities. c PARAMETER (MAXDIM=100,MAXTR=4*MAXDIM) c COMMON /FRBLK/WK COMMON /LNUMS/IPT INTEGER*4 IPT(MAXTR),IIPT,L,I,J,IR REAL*8 WK(MAXTR),WWK C L=N/2+1 IR=N 10 CONTINUE IF(L.GT.NSTART)THEN L=L-1 WWK=WK(L) IIPT=IPT(L) ELSE WWK=WK(IR) IIPT=IPT(IR) WK(IR)=WK(1) IPT(IR)=IPT(1) IR=IR-1 IF(IR.EQ.NSTART)THEN WK(1)=WWK IPT(1)=IIPT RETURN ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(WK(J).LT.WK(J+1))J=J+1 ENDIF IF(WWK.LT.WK(J))THEN WK(I)=WK(J) IPT(I)=IPT(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF WK(I)=WWK IPT(I)=IIPT GO TO 10 c RETURN END C C------------------------------------------------------------------------------ C------------------------------------------------------------------------------