1 SUBROUTINE FREQCY(FMATRX,FREQ,CNORML,REDMAS,TRAVEL,EORC,DELDIP) 2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 3 INCLUDE 'SIZES' 4 DIMENSION FMATRX(*), REDMAS(*), FREQ(*), CNORML(*), TRAVEL(*) 5 +,DELDIP(3,*) 6 LOGICAL EORC 7********************************************************************* 8* 9* FRCE CALCULATES THE FORCE CONSTANTS AND VIBRATIONAL FREQUENCIES 10* FOR A MOLECULE. IT USES THE ISOTOPIC MASSES TO WEIGHT THE 11* FORCE MATRIX 12* 13* ON INPUT FMATRX = FORCE MATRIX, OF SIZE NUMAT*3*(NUMAT*3+1)/2. 14* 15********************************************************************* 16 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 17 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 18 2 NCLOSE,NOPEN,NDUMY,FRACT 19 COMMON /NLLCOM/ FMAT2D(2*MAXPAR**2), VEC(MAXPAR**2) 20 COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120), NENT 21 COMMON /ATMASS/ ATMASS(NUMATM) 22 COMMON /KEYWRD/ KEYWRD 23 COMMON /SCRACH/ OLDF(MAXPAR**2) 24 COMMON /WORK1 / DUMMY1(NPULAY*4), DUMMY2(NPULAY*2), 25 . DUMMY3(NPULAY*2), ALBAND(NPULAY*13) 26 27 CHARACTER KEYWRD*241 28 DIMENSION WTMASS(MAXPAR), SHIFT(6), SEC(MAXPAR**2) 29 COMPLEX SEC, VEC 30 EQUIVALENCE (SEC,OLDF) 31 SAVE FACT 32 DATA FACT/6.023D23/ 33C 34C CONVERSION FACTOR FOR SPEED OF LIGHT AND 2 PI. 35C 36 C2PI=1.D0/(2.998D10*3.141592653598D0*2.D0) 37C NOW TO CALCULATE THE VIBRATIONAL FREQUENCIES 38C 39C FIND CONVERSION CONSTANTS FOR MASS WEIGHTED SYSTEM 40 IF(INDEX(KEYWRD,' GROUP').NE.0) THEN 41 CALL SYMT(FMATRX, DELDIP) 42 IF(INDEX(KEYWRD,' FREQCY').NE.0)THEN 43 WRITE(6,'(A)')' SYMMETRIZED HESSIAN MATRIX' 44C# I=-N3 45C# CALL VECPRT(FMATRX,I) 46 47C 48C THE FORCE MATRIX IS PRINTED AS AN ATOM-ATOM MATRIX RATHER THAN 49C AS A 3N*3N MATRIX, AS THE 3N MATRIX IS VERY CONFUSING! 50C 51 IJ=0 52 IU=0 53 DO 159 I=1,NUMAT 54 IL=IU+1 55 IU=IL+2 56 IM1=I-1 57 JU=0 58 DO 149 J=1,IM1 59 JL=JU+1 60 JU=JL+2 61 SUM=0.D0 62C$DOIT ASIS 63 DO 139 II=IL,IU 64C$DOIT ASIS 65 DO 139 JJ=JL,JU 66 139 SUM=SUM+FMATRX((II*(II-1))/2+JJ)**2 67 IJ=IJ+1 68 149 CNORML(IJ)=SQRT(SUM) 69 IJ=IJ+1 70 159 CNORML(IJ)=SQRT( 71 1FMATRX(((IL+0)*(IL+1))/2)**2+ 72 2FMATRX(((IL+1)*(IL+2))/2)**2+ 73 3FMATRX(((IL+2)*(IL+3))/2)**2+2.D0*( 74 4FMATRX(((IL+1)*(IL+2))/2-1)**2+ 75 5FMATRX(((IL+2)*(IL+3))/2-2)**2+ 76 6FMATRX(((IL+2)*(IL+3))/2-1)**2)) 77 I=-NUMAT 78 CALL VECPRT(CNORML,I) 79 ENDIF 80 ENDIF 81 N3=3*NUMAT 82 L=0 83 DO 10 I=1,NUMAT 84 WEIGHT=1.4142136D0/SQRT(ATMASS(I)) 85 WTMASS(L+1)=WEIGHT 86 WTMASS(L+2)=WEIGHT 87 WTMASS(L+3)=WEIGHT 88 L=L+3 89 10 WTMASS(L)=WEIGHT 90C CONVERT TO MASS WEIGHTED FMATRX 91 LINEAR=0 92 DO 20 I=1,N3 93 DO 20 J=1,I 94 LINEAR=LINEAR+1 95 OLDF(LINEAR)= FMATRX(LINEAR)*1.D5 96 20 FMATRX(LINEAR)=FMATRX(LINEAR)*WTMASS(I)*WTMASS(J) 97C 98C 1.D5 IS TO CONVERT FROM MILLIDYNES/ANGSTROM TO DYNES/CM. 99C 100C DIAGONALIZE 101 I=INDEX(KEYWRD,' K=') 102 IF(I.NE.0)THEN 103C 104C GO INTO BRILLOUIN ZONE MODE 105C 106 STEP=READA(KEYWRD,I) 107 MONO3=READA(KEYWRD(I:),INDEX(KEYWRD(I:),','))*3 108 CALL BRLZON(FMATRX, FMAT2D, N3, SEC, VEC, ALBAND, MONO3, STEP,1 109 1) 110 RETURN 111 ENDIF 112 CALL FRAME(FMATRX,NUMAT,1, SHIFT) 113 CALL RSP(FMATRX,N3,N3,FREQ,CNORML) 114 DO 30 I=1,N3 115 J=(FREQ(I)+50.D0)*0.01D0 116 30 FREQ(I)=FREQ(I)-DBLE(J*100) 117 DO 40 I=1,N3 118 40 FREQ(I)=FREQ(I)*1.D5 119C 120C CALCULATE REDUCED MASSES, STORE IN REDMAS 121C 122 DO 80 I=1,N3 123 II=(I-1)*N3 124 SUM=0.D0 125 DO 70 J=1,N3 126 JII=J+II 127 JJ=(J*(J-1))/2 128 DO 50 K=1,J 129 50 SUM=SUM+CNORML(JII)*OLDF(JJ+K)*CNORML(K+II) 130 DO 60 K=J+1,N3 131 60 SUM=SUM+CNORML(JII)*OLDF((K*(K-1))/2+J)*CNORML(K+II) 132 70 CONTINUE 133 SUM1=SUM*2.D0 134 IF(ABS(FREQ(I)).GT.ABS(SUM)*1.D-20) THEN 135 SUM=1.D0*SUM/FREQ(I) 136 ELSE 137 SUM=0.D0 138 ENDIF 139 FREQ(I)=SIGN(SQRT(FACT*ABS(FREQ(I)))*C2PI,FREQ(I)) 140 IF(ABS(FREQ(I)).LT.ABS(SUM1)*1.D+20) THEN 141 SUM1=SQRT(ABS(FREQ(I)/(SUM1*1.D-5))) 142 ELSE 143 SUM1=0.D0 144 ENDIF 145 IF(SUM.LT.0.D0.OR.SUM.GT.100)SUM=0.D0 146C 147C 0.0063024=SQRT(2*A*B*C/N) WHERE 148C A=1.196D8 = CONVERSION OF CM**(-1) TO (ERGS = DYNE.ANGSTROMS) 149C B=1000.0 = MILLIDYNES TO DYNES 150C C=1.D8 = CENTIMETERS TO ANGSTROMS 151C N=6.02205D23 = AVOGADRO'S NUMBER 152 TRAVEL(I)=SUM1*0.0063024D0 153 IF(TRAVEL(I).GT.1.D0)TRAVEL(I)=0.D0 154C# WRITE(6,*)TRAVEL(I) 155 80 REDMAS(I)=SUM 156 IF(INDEX(KEYWRD,' GROUP').NE.0) CALL SYMA(FREQ, CNORML) 157 IF(EORC) THEN 158C 159C CONVERT NORMAL VECTORS TO CARTESIAN COORDINATES 160C (DELETED) AND NORMALIZE SO THAT THE TOTAL MOVEMENT IS 1.0 ANGSTROM. 161C 162 IJ=0 163 DO 110 I=1,N3 164 SUM=0.D0 165 J=0 166 DO 90 JJ=1,NUMAT 167 SUM1=0.D0 168 CNORML(IJ+1)=CNORML(IJ+1)*WTMASS(J+1) 169 SUM1=SUM1+CNORML(IJ+1)**2 170C 171 CNORML(IJ+2)=CNORML(IJ+2)*WTMASS(J+2) 172 SUM1=SUM1+CNORML(IJ+2)**2 173C 174 CNORML(IJ+3)=CNORML(IJ+3)*WTMASS(J+3) 175 SUM1=SUM1+CNORML(IJ+3)**2 176C 177 J=J+3 178 IJ=IJ+3 179 90 SUM=SUM+SQRT(SUM1) 180 SUM=1.D0/SQRT(SUM) 181 IJ=IJ-N3 182 DO 100 J=1,N3 183 IJ=IJ+1 184 100 CNORML(IJ)=CNORML(IJ)*SUM 185 110 CONTINUE 186C 187C RETURN HESSIAN IN MILLIDYNES/ANGSTROM IN FMATRX 188C 189 DO 120 I=1,LINEAR 190 120 FMATRX(I)=OLDF(I)*1.D-5 191 ELSE 192C 193C RETURN HESSIAN AS MASS-WEIGHTED FMATRIX 194 LINEAR=0 195C 196 DO 130 I=1,N3 197 DO 130 J=1,I 198 LINEAR=LINEAR+1 199 130 FMATRX(LINEAR)=OLDF(LINEAR)*1.D-5*WTMASS(I)*WTMASS(J) 200 ENDIF 201 RETURN 202 END 203