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