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