SUBROUTINE ANAVIB(COORD,EIGS,DIPT,N3,VIBS,RIJ,HESS,TRAVEL,REDMAS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE 'SIZES' DIMENSION COORD(3,NUMATM),EIGS(N3),VIBS(N3,N3), DIPT(N3), 1RIJ(MAXHES), TRAVEL(N3), HESS(*), REDMAS(*) COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 2 NCLOSE,NOPEN,NDUMY,FRACT COMMON /ELEMTS/ ELEMNT(107) LOGICAL VIB1, VIB2, VIB3, VIB4, VIB5 COMMON /FOKMAT/ F(MPACK*2) CHARACTER*2 ELEMNT DIMENSION VANRAD(107), IJF(10), FIJ(10) SAVE VANRAD, FACT DATA VANRAD/ 1 0.32,0.93, 2 1.23, 0.90, 0.82, 0.77, 0.75, 0.73, 0.72, 0.71, 3 1.54, 1.36, 1.18, 1.11, 1.06, 1.02, 0.99, 0.98, 4 2.03, 1.74, 1.44, 1.32, 1.22, 1.18, 1.17, 1.17, 1.16, 5 1.15, 1.17, 1.25, 1.26, 1.22, 1.20, 1.16, 1.14, 1.12, 6 2.16, 1.91, 1.62, 1.45, 1.34, 1.30, 1.27, 1.25, 1.25, 7 1.28, 1.34, 1.48, 1.44, 1.41, 1.40, 1.36, 1.33, 1.31, 8 2.35, 1.98, 1.69, 9 1.65, 1.65, 1.64, 1.63, 1.62, 1.85, 1.61, 1.59, 1.59, 1.58, 1 1.57, 1.56, 1.56, 1.56, 2 1.44, 1.34, 1.30, 1.28, 1.26, 1.27, 1.30, 1.34, 3 1.49, 1.48, 1.47, 1.46, 1.46, 1.45,1.45,21*1.45/ DATA FACT/6.023D23/ C2PI=1.D0/(2.998D10*3.141592653598D0*2.D0) N3=NUMAT*3 C C COMPUTE INTERATOMIC DISTANCES. C L=0 DO 10 I=1,NUMAT DO 10 J=1,I-1 L=L+1 10 RIJ(L)=SQRT((COORD(1,J)-COORD(1,I))**2+ 1 (COORD(2,J)-COORD(2,I))**2+ 2 (COORD(3,J)-COORD(3,I))**2)+1.D-10 C C ANALYSE VIBRATIONS C WRITE(6,'(''1'',//10X,''DESCRIPTION OF VIBRATIONS'',/)') ILINE=0 DO 120 K=1,N3 IF(ABS(EIGS(K)).LT.50) GOTO 120 VIB1=.TRUE. VIB2=.TRUE. VIB3=.TRUE. VIB4=.TRUE. VIB5=.TRUE. J3=0 L=0 TOT=0.D0 LINEAR=0 J1=-2 DO 70 J=1,NUMAT J1=J1+3 I1=-2 DO 70 I=1,J-1 I1=I1+3 VDW=(VANRAD(NAT(I))+VANRAD(NAT(J)))*1.5 L=L+1 F(L)=0.D0 IF( RIJ(L) .LT. VDW) THEN C C CALCULATE ENERGY TERM BETWEEN THE TWO ATOMS C EAB=0.D0 DO 20 JJ=J1,J1+2 DO 20 II=I1,I1+2 20 EAB=EAB+VIBS(JJ,K)*HESS((JJ*(JJ-1))/2+II)*VIBS(II,K) EB=0.D0 DO 40 JJ=J1,J1+2 DO 30 II=J1,JJ 30 EB=EB+VIBS(JJ,K)*HESS((JJ*(JJ-1))/2+II)*VIBS(II,K)* 12.D0 40 EB=EB-VIBS(JJ,K)*HESS((JJ*(JJ+1))/2)*VIBS(JJ,K) EA=0.D0 DO 60 JJ=I1,I1+2 DO 50 II=I1,JJ 50 EA=EA+VIBS(JJ,K)*HESS((JJ*(JJ-1))/2+II)*VIBS(II,K)* 12.D0 60 EA=EA-VIBS(JJ,K)*HESS((JJ*(JJ+1))/2)*VIBS(JJ,K) LINEAR=LINEAR+1 F(L)=EA+EAB*2.D0+EB TOT=TOT+F(L) ENDIF 70 CONTINUE C C NOW TO SORT F INTO DECENDING ORDER C DO 90 I=1,10 SUM=-100.D0 DO 80 J=1,L IF(ABS(F(J)).GT.SUM)THEN JJ=J SUM=ABS(F(J)) ENDIF 80 CONTINUE IF(SUM.LT.0.D0)GOTO 100 FIJ(I)=SUM F(JJ)=-1.D-5 IJF(I)=JJ C# WRITE(6,*)FIJ(I),IJF(I) 90 CONTINUE I=10 100 LINEAR=I SUM=1.D0/(TOT+1.D-8) DO 110 IJ=1,LINEAR J=0.5D0*(0.99D0+SQRT(1.D0+8.D0*IJF(IJ))) I=IJF(IJ)-(J*(J-1))/2 J=J+1 XJ=COORD(1,J) YJ=COORD(2,J) ZJ=COORD(3,J) J1=3*J-2 J2=J1+1 J3=J2+1 I3=0 XI=COORD(1,I) YI=COORD(2,I) ZI=COORD(3,I) I1=3*I-2 I2=I1+1 I3=I2+1 X= VIBS(J1,K)-VIBS(I1,K) Y= VIBS(J2,K)-VIBS(I2,K) Z= VIBS(J3,K)-VIBS(I3,K) E=FIJ(IJ)*SUM*100.D0 C# IF(ABS(E).GT.110)GOTO 120 SHIFT=X*X+Y*Y+Z*Z+1.D-30 IF(ABS(E) .GT. 10.D0.OR.IJ.LT.5.AND.ABS(E).GT.0.1) THEN SHIFT=SQRT(SHIFT) RADIAL=((X*(XI-XJ)+Y*(YI-YJ)+Z*(ZI-ZJ)) 1 /(SHIFT*RIJ(IJF(IJ))))**2*100.D0 IF (VIB1) THEN WRITE(6,'(/,'' VIBRATION'',I4,'' ATOM PAIR 1 '','' ENERGY CONTRIBUTION RADIAL'')')K ANS=100.D0*SQRT(FIJ(IJ)*1.D5*6.023D23)/(2.998D10*3.141 159D0*2.D0)/EIGS(K) ANS=MIN(999.9D0,MAX(-99.9D0,ANS)) WRITE(6,'('' FREQ. '',F9.2,6X,1A2,I2, 1'' -- '',A2,I2, '' '',F6.1,''% ('',F5.1,''%)'',F18.1 2,''%'')') 3EIGS(K),ELEMNT(NAT(I)),I,ELEMNT(NAT(J)),J,E,ANS,RADIAL C# WRITE(6,*)ANS VIB1=.FALSE. ELSEIF (VIB2) THEN VIB2=.FALSE. WRITE(6,'('' T-DIPOLE'',F9.4,6X,1A2,I2, 1'' -- '',A2,I2, '' '',F6.1,''%'',F27.1,''%'')') 2DIPT(K),ELEMNT(NAT(I)),I,ELEMNT(NAT(J)),J,E,RADIAL ELSEIF (VIB3) THEN VIB3=.FALSE. WRITE(6,'('' TRAVEL '',F9.4,6X,1A2,I2, 1'' -- '',A2,I2, '' '',F6.1,''%'',F27.1,''%'')') 2TRAVEL(K),ELEMNT(NAT(I)),I,ELEMNT(NAT(J)),J,E,RADIAL ELSEIF (VIB4) THEN VIB4=.FALSE. WRITE(6,'('' RED. MASS'',F8.4,6X,1A2,I2, 1'' -- '',A2,I2, '' '',F6.1,''%'',F27.1,''%'')') 2REDMAS(K),ELEMNT(NAT(I)),I,ELEMNT(NAT(J)),J,E,RADIAL ELSE ILINE=ILINE+1 WRITE(6,'('' '',1A2,I2, 1'' -- '',A2,I2, '' '',F6.1,''%'',F27.1,''%'')') 2ELEMNT(NAT(I)),I,ELEMNT(NAT(J)),J,E,RADIAL ENDIF ENDIF 110 CONTINUE ILINE=ILINE+6 IF(VIB1)WRITE(6,'(/,'' VIBRATION'',I4)')K IF(VIB1)WRITE(6,'( '' FREQ. '',F8.2)')EIGS(K) IF(VIB2)WRITE(6,'( '' T-DIPOLE '',F8.4)')DIPT(K) IF(VIB3)WRITE(6,'( '' TRAVEL '',F8.4)')TRAVEL(K) IF(VIB4)WRITE(6,'( '' RED. MASS'',F8.4)')REDMAS(K) IF(ILINE.GT.52)THEN ILINE=0 WRITE(6,'(''1'')') ENDIF 120 CONTINUE RETURN END