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
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
|
SUBROUTINE REACT1(ESCF)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
REAL PASTOR, PBSTOR
INCLUDE 'SIZES'
COMMON /GEOM / GEO(3,NUMATM), XCOORD(3,NUMATM)
DIMENSION GEOA(3,NUMATM), GEOVEC(3,NUMATM),
1 PASTOR(MPACK),
2 PBSTOR(MPACK), XOLD(MAXPAR), GROLD(MAXPAR)
COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
1 NA(NUMATM), NB(NUMATM), NC(NUMATM)
COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
COMMON /GEOSYM/ NDEP,LOCPAR(MAXPAR),IDEPFN(MAXPAR),LOCDEP(MAXPAR)
COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
COMMON /GRADNT/ GRAD(MAXPAR),GNORM
COMMON /ISTOPE/ AMS(107)
COMMON /GRAVEC/ COSINE
COMMON /KEYWRD/ KEYWRD
COMMON /MESAGE/ IFLEPO,ISCF
1 /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
2 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
3 NCLOSE,NOPEN,NDUMY,FRACT
COMMON /REACTN/ STEP, GEOA, GEOVEC,CALCST
LOGICAL GRADNT, FINISH, XYZ, INT, GOK(2)
SAVE GRADNT, FINISH, XYZ,INT, GOK
************************************************************************
*
* REACT1 DETERMINES THE TRANSITION STATE OF A CHEMICAL REACTION.
*
* REACT WORKS BY USING TWO SYSTEMS SIMULTANEOUSLY, THE HEATS OF
* FORMATION OF BOTH ARE CALCULATED, THEN THE MORE STABLE ONE
* IS MOVED IN THE DIRECTION OF THE OTHER. AFTER A STEP THE
* ENERGIES ARE COMPARED, AND THE NOW LOWER-ENERGY FORM IS MOVED
* IN THE DIRECTION OF THE HIGHER-ENERGY FORM. THIS IS REPEATED
* UNTIL THE SADDLE POINT IS REACHED.
*
* IF ONE FORM IS MOVED 3 TIMES IN SUCCESSION, THEN THE HIGHER ENERGY
* FORM IS RE-OPTIMIZED WITHOUT SHORTENING THE DISTANCE BETWEEN THE TWO
* FORMS. THIS REDUCES THE CHANCE OF BEING CAUGHT ON THE SIDE OF A
* TRANSITION STATE.
*
************************************************************************
DIMENSION IDUM1(NUMATM), IDUM2(NUMATM), XSTORE(MAXPAR),
1IDUM3(NUMATM), COORD(3,NUMATM), IROT(2,3)
DIMENSION IDUMMY(3*NUMATM)
SAVE IROT
CHARACTER*241 KEYWRD
EQUIVALENCE (IDUMMY,COORD)
DATA IROT/1,2,1,3,2,3/
GOLD=0.D0
LINEAR=0
IFLAG=1
GOK(1)=.FALSE.
GOK(2)=.FALSE.
XYZ=(INDEX(KEYWRD,' XYZ') .NE. 0)
GRADNT=(INDEX(KEYWRD,'GRAD') .NE. 0)
I=(INDEX(KEYWRD,' BAR'))
STEPMX=0.15D0
IF(I.NE.0) STEPMX=READA(KEYWRD,I)
MAXSTP=1000
C
C READ IN THE SECOND GEOMETRY.
C
IF(XYZ) THEN
CALL GETGEO(5,LABELS,GEOA,LOC,NA,NB,NC,AMS,NATOMS,INT)
ELSE
CALL GETGEO(5,IDUM1,GEOA,IDUMMY,
1 IDUM1,IDUM2,IDUM3,AMS,NATOMS,INT)
C
C IF INTERNAL COORDINATES ARE TO BE USED, CHECK THE CONNECTIVITY
C
L=0
DO 10 I=1,NATOMS
IF(IDUM1(I).NE.NA(I))THEN
L=L+1
IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY
1'')')
WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE BOND LAB
1ELS ARE DIFFERENT: ',IDUM1(I),' AND',NA(I)
ENDIF
IF(IDUM2(I).NE.NB(I))THEN
L=L+1
IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY
1'')')
WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE BOND ANG
1LE LABELS ARE DIFFERENT:',IDUM2(I),' AND',NB(I)
ENDIF
IF(IDUM3(I).NE.NC(I))THEN
L=L+1
IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY
1'')')
WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE DIHEDRAL
1 LABELS ARE DIFFERENT: ',IDUM3(I),' AND',NC(I)
ENDIF
10 CONTINUE
IF(L.NE.0)WRITE(6,'(10X,A)')' CORRECT BEFORE RESUBMISSION'
IF(L.NE.0)STOP
ENDIF
TIME0= SECOND()
C
C SWAP FIRST AND SECOND GEOMETRIES AROUND
C SO THAT GEOUT CAN OUTPUT DATA ON SECOND GEOMETRY.
C
NUMAT2=0
DO 20 I=1,NATOMS
IF(LABELS(I).NE.99) NUMAT2=NUMAT2+1
X=GEOA(1,I)
GEOA(1,I)=GEO(1,I)
GEO(1,I)=X
X=GEOA(2,I)*0.0174532925D0
GEOA(2,I)=GEO(2,I)
GEO(2,I)=X
X=GEOA(3,I)*0.0174532925D0
GEOA(3,I)=GEO(3,I)
GEO(3,I)=X
20 CONTINUE
IF(NUMAT2.NE.NUMAT) THEN
WRITE(6,'(//10X,'' NUMBER OF ATOMS IN SECOND SYSTEM IS '',
1''INCORRECT'',/)')
WRITE(6,'('' NUMBER OF ATOMS IN FIRST SYSTEM ='',I4)')NUMAT
WRITE(6,'('' NUMBER OF ATOMS IN SECOND SYSTEM ='',I4)')NUMAT2
GOTO 280
ENDIF
WRITE(6,'(//10X,'' GEOMETRY OF SECOND SYSTEM'',/)')
IF(NDEP.NE.0) CALL SYMTRY
CALL GEOUT(1)
C
C CONVERT TO CARTESIAN, IF NECESSARY
C
IF( XYZ )THEN
CALL GMETRY(GEO,COORD)
SUMX=0.D0
SUMY=0.D0
SUMZ=0.D0
DO 30 J=1,NUMAT
SUMX=SUMX+COORD(1,J)
SUMY=SUMY+COORD(2,J)
30 SUMZ=SUMZ+COORD(3,J)
SUMX=SUMX/NUMAT
SUMY=SUMY/NUMAT
SUMZ=SUMZ/NUMAT
DO 40 J=1,NUMAT
GEO(1,J)=COORD(1,J)-SUMX
GEO(2,J)=COORD(2,J)-SUMY
40 GEO(3,J)=COORD(3,J)-SUMZ
WRITE(6,'(//,'' CARTESIAN GEOMETRY OF FIRST SYSTEM'',//)')
WRITE(6,'(3F14.5)')((GEO(J,I),J=1,3),I=1,NUMAT)
SUMX=0.D0
SUMY=0.D0
SUMZ=0.D0
DO 50 J=1,NUMAT
SUMX=SUMX+GEOA(1,J)
SUMY=SUMY+GEOA(2,J)
50 SUMZ=SUMZ+GEOA(3,J)
SUM=0.D0
SUMX=SUMX/NUMAT
SUMY=SUMY/NUMAT
SUMZ=SUMZ/NUMAT
DO 60 J=1,NUMAT
GEOA(1,J)=GEOA(1,J)-SUMX
GEOA(2,J)=GEOA(2,J)-SUMY
GEOA(3,J)=GEOA(3,J)-SUMZ
SUM=SUM+(GEO(1,J)-GEOA(1,J))**2
1 +(GEO(2,J)-GEOA(2,J))**2
2 +(GEO(3,J)-GEOA(3,J))**2
60 CONTINUE
DO 110 L=3,1,-1
C
C DOCKING IS DONE IN STEPS OF 16, 4, AND 1 DEGREES AT A TIME.
C
CA=COS(4.D0**(L-1)*0.01745329D0)
SA=SQRT(ABS(1.D0-CA**2))
DO 100 J=1,3
IR=IROT(1,J)
JR=IROT(2,J)
DO 90 I=1,10
SUMM=0.D0
DO 70 K=1,NUMAT
X = CA*GEOA(IR,K)+SA*GEOA(JR,K)
GEOA(JR,K)=-SA*GEOA(IR,K)+CA*GEOA(JR,K)
GEOA(IR,K)=X
SUMM=SUMM+(GEO(1,K)-GEOA(1,K))**2
1 +(GEO(2,K)-GEOA(2,K))**2
2 +(GEO(3,K)-GEOA(3,K))**2
70 CONTINUE
IF(SUMM.GT.SUM) THEN
IF(I.GT.1)THEN
SA=-SA
DO 80 K=1,NUMAT
X = CA*GEOA(IR,K)+SA*GEOA(JR,K)
GEOA(JR,K)=-SA*GEOA(IR,K)+CA*GEOA(JR,K)
GEOA(IR,K)=X
80 CONTINUE
GOTO 100
ENDIF
SA=-SA
ENDIF
90 SUM=SUMM
100 CONTINUE
110 CONTINUE
WRITE(6,'(//,'' CARTESIAN GEOMETRY OF SECOND SYSTEM'',//)')
WRITE(6,'(3F14.5)')((GEOA(J,I),J=1,3),I=1,NUMAT)
WRITE(6,'(//,'' "DISTANCE":'',F13.6)')SUM
WRITE(6,'(//,'' REACTION COORDINATE VECTOR'',//)')
WRITE(6,'(3F14.5)')((GEOA(J,I)-GEO(J,I),J=1,3),I=1,NUMAT)
NA(1)=99
J=0
NVAR=0
DO 130 I=1,NATOMS
IF(LABELS(I).NE.99)THEN
J=J+1
DO 120 K=1,3
NVAR=NVAR+1
LOC(2,NVAR)=K
120 LOC(1,NVAR)=J
LABELS(J)=LABELS(I)
ENDIF
130 CONTINUE
NATOMS=NUMAT
ENDIF
C
C XPARAM HOLDS THE VARIABLE PARAMETERS FOR GEOMETRY IN GEO
C XOLD HOLDS THE VARIABLE PARAMETERS FOR GEOMETRY IN GEOA
C
IF(NVAR.EQ.0)THEN
WRITE(6,'(///10X,''THERE ARE NO VARIABLES IN THE SADDLE'',
1'' CALCULATION!'')')
STOP
ENDIF
SUM=0.D0
DO 140 I=1,NVAR
GROLD(I)=1.D0
XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
XOLD(I)=GEOA(LOC(2,I),LOC(1,I))
140 SUM=SUM+(XPARAM(I)-XOLD(I))**2
STEP0=SQRT(SUM)
IF(STEP0.LT.1.D-5)THEN
WRITE(6,'(//,3(5X,A,/))')' BOTH GEOMETRIES ARE IDENTICAL',
1' A SADDLE CALCULATION INVOLVES A REACTANT AND A PRODUCT',
2' THESE MUST BE DIFFERENT GEOMETRIES'
STOP
ENDIF
ONE=1.D0
DELL=0.1D0
EOLD=-2000.D0
TIME1=SECOND()
SWAP=0
DO 240 ILOOP=1,MAXSTP
WRITE(6,'('' '',40(''*+''))')
C
C THIS METHOD OF CALCULATING 'STEP' IS QUITE ARBITARY, AND NEEDS
C TO BE IMPROVED BY INTELLIGENT GUESSWORK!
C
IF (GNORM.LT.1.D-3)GNORM=1.D-3
STEP=MIN(SWAP,0.5D0, 6.D0/GNORM, DELL,STEPMX*STEP0+0.005D0)
STEP=MIN(0.2D0,STEP/STEP0)*STEP0
SWAP=SWAP+1.D0
DELL=DELL+0.1D0
WRITE(6,'('' BAR SHORTENED BY'',F12.7,'' PERCENT'')')
1STEP/STEP0*100.D0
STEP0=STEP0-STEP
IF(STEP0.LT.0.01D0) GOTO 250
STEP=STEP0
DO 150 I=1,NVAR
150 XSTORE(I)=XPARAM(I)
CALL FLEPO(XPARAM, NVAR, ESCF)
IF(LINEAR.EQ.0)THEN
LINEAR=(NORBS*(NORBS+1))/2
DO 160 I=1,LINEAR
PASTOR(I)=PA(I)
160 PBSTOR(I)=PB(I)
ENDIF
DO 170 I=1,NVAR
170 XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
IF(IFLAG.EQ.1)THEN
WRITE(6,'(//10X,''FOR POINT'',I3,'' SECOND STRUCTURE'')')ILO
1OP
ELSE
WRITE(6,'(//10X,''FOR POINT'',I3,'' FIRST STRUCTURE'')')ILO
1OP
ENDIF
WRITE(6,'('' DISTANCE A - B '',F12.6)')STEP
C
C NOW TO CALCULATE THE "CORRECT" GRADIENTS, SWITCH OFF 'STEP'.
C
STEP=0.D0
DO 180 I=1,NVAR
180 GRAD(I)=GROLD(I)
CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.)
DO 190 I=1,NVAR
190 GROLD(I)=GRAD(I)
IF (GRADNT) THEN
WRITE(6,'('' ACTUAL GRADIENTS OF THIS POINT'')')
WRITE(6,'(8F10.4)')(GRAD(I),I=1,NVAR)
ENDIF
WRITE(6,'('' HEAT '',F12.6)')FUNCT1
GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
WRITE(6,'('' GRADIENT NORM '',F12.6)')GNORM
COSINE=COSINE*ONE
WRITE(6,'('' DIRECTION COSINE'',F12.6)')COSINE
CALL GEOUT(6)
IF(SWAP.GT.2.9D0 .OR. ILOOP .GT. 3 .AND. COSINE .LT. 0.D0
1 .OR. ESCF .GT. EOLD)
2 THEN
IF(SWAP.GT.2.9D0)THEN
SWAP=0.D0
ELSE
SWAP=0.5D0
ENDIF
C
C SWAP REACTANT AND PRODUCT AROUND
C
FINISH=(GOK(1).AND.GOK(2) .AND. COSINE .LT. 0.D0)
IF(FINISH) THEN
WRITE(6,'(//10X,'' BOTH SYSTEMS ARE ON THE SAME SIDE OF T
1HE '',''TRANSITION STATE -'',/10X,'' GEOMETRIES OF THE SYSTEMS'',
2'' ON EACH SIDE OF THE T.S. ARE AS FOLLOWS'')')
DO 200 I=1,NVAR
200 XPARAM(I)=XSTORE(I)
CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.)
WRITE(6,'(//10X,'' GEOMETRY ON ONE SIDE OF THE TRANSITION
1'','' STATE'')')
CALL WRITMO(TIME0,FUNCT1)
ENDIF
TIME2=SECOND()
WRITE(6,'('' TIME='',F9.2)')TIME2-TIME1
TIME1=TIME2
WRITE(6,'('' REACTANTS AND PRODUCTS SWAPPED AROUND'')')
IFLAG=1-IFLAG
ONE=-1.D0
EOLD=ESCF
SUM=GOLD
GOLD=GNORM
I=1.7D0+ONE*0.5D0
IF(GNORM.GT.10.D0)GOK(I)=.TRUE.
GNORM=SUM
DO 210 I=1,NATOMS
DO 210 J=1,3
X=GEO(J,I)
GEO(J,I)=GEOA(J,I)
210 GEOA(J,I)=X
DO 220 I=1,NVAR
X=XOLD(I)
XOLD(I)=XPARAM(I)
220 XPARAM(I)=X
C
C
C SWAP AROUND THE DENSITY MATRICES.
C
DO 230 I=1,LINEAR
X=PASTOR(I)
PASTOR(I)=PA(I)
PA(I)=X
X=PBSTOR(I)
PBSTOR(I)=PB(I)
PB(I)=X
P(I)=PA(I)+PB(I)
230 CONTINUE
IF(FINISH) GOTO 250
ELSE
ONE=1.D0
ENDIF
240 CONTINUE
250 CONTINUE
WRITE(6,'('' AT END OF REACTION'')')
GOLD=SQRT(DOT(GRAD,GRAD,NVAR))
CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.)
GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
DO 260 I=1,NVAR
260 GROLD(I)=XPARAM(I)
CALL WRITMO(TIME0,FUNCT1)
*
* THE GEOMETRIES HAVE (A) BEEN OPTIMIZED CORRECTLY, OR
* (B) BOTH ENDED UP ON THE SAME SIDE OF THE T.S.
*
* TRANSITION STATE LIES BETWEEN THE TWO GEOMETRIES
*
C1=GOLD/(GOLD+GNORM)
C2=1.D0-C1
WRITE(6,'('' BEST ESTIMATE GEOMETRY OF THE TRANSITION STATE'')')
WRITE(6,'(//10X,'' C1='',F8.3,''C2='',F8.3)')C1,C2
DO 270 I=1,NVAR
270 XPARAM(I)=C1*GROLD(I)+C2*XOLD(I)
STEP=0.D0
CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.)
CALL WRITMO(TIME0,FUNCT1)
280 RETURN
END
|
