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
388
389
390
391
392
|
SUBROUTINE GETGEO(IREAD,LABELS,GEO,LOPT,NA,NB,NC,AMS,NATOMS,INT)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'SIZES'
DIMENSION GEO(3,*),NA(*),NB(*),NC(*),AMS(*), LOPT(3,*)
1,LABELS(*)
LOGICAL INT
************************************************************************
*
* GETGEO READS IN THE GEOMETRY. THE ELEMENT IS SPECIFIED BY IT'S
* CHEMICAL SYMBOL, OR, OPTIONALLY, BY IT'S ATOMIC NUMBER.
*
* ON INPUT IREAD = CHANNEL NUMBER FOR READ, NORMALLY 5
* AMS = DEFAULT ATOMIC MASSES.
*
* ON OUTPUT LABELS = ATOMIC NUMBERS OF ALL ATOMS, INCLUDING DUMMIES.
* GEO = INTERNAL COORDINATES, IN ANGSTROMS, AND DEGREES.
* LOPT = INTEGER ARRAY, A '1' MEANS OPTIMIZE THIS PARAMETER,
* '0' MEANS DO NOT OPTIMIZE, AND A '-1' LABELS THE
* REACTION COORDINATE.
* NA = INTEGER ARRAY OF ATOMS (SEE DATA INPUT)
* NB = INTEGER ARRAY OF ATOMS (SEE DATA INPUT)
* NC = INTEGER ARRAY OF ATOMS (SEE DATA INPUT)
* ATMASS = ATOMIC MASSES OF ATOMS.
************************************************************************
COMMON /PATH / IDUM(2),REACT(3,66), DUMM1,DUMM2
COMMON /SIMBOL/ SIMBOL(MAXPAR)
COMMON /ATMASS/ ATMASS(NUMATM)
COMMON /ATOMTX/ LTXT, TXTATM(NUMATM)
COMMON /KEYWRD/ KEYWRD
DIMENSION ISTART(40), XYZ(3,NUMATM), VALUE(40)
LOGICAL LEADSP, IRCDRC
CHARACTER KEYWRD*241, TXTATM*8, SIMBOL*10, LTXT*1
CHARACTER ELEMNT(107)*2, LINE*80, SPACE*1, NINE*1,ZERO*1,
1TAB*1, COMMA*1, STRING*80, ELE*2, TURN*1
SAVE ELEMNT, COMMA, SPACE, NINE, ZERO
DATA (ELEMNT(I),I=1,107)/'H','HE',
1 'LI','BE','B','C','N','O','F','NE',
2 'NA','MG','AL','SI','P','S','CL','AR',
3 'K','CA','SC','TI','V','CR','MN','FE','CO','NI','CU',
4 'ZN','GA','GE','AS','SE','BR','KR',
5 'RB','SR','Y','ZR','NB','MO','TC','RU','RH','PD','AG',
6 'CD','IN','SN','SB','TE','I','XE',
7 'CS','BA','LA','CE','PR','ND','PM','SM','EU','GD','TB','DY',
8 'HO','ER','TM','YB','LU','HF','TA','W','RE','OS','IR','PT',
9 'AU','HG','TL','PB','BI','PO','AT','RN',
1 'FR','RA','AC','TH','PA','U','NP','PU','AM','CM','BK','CF','XX',
2 'FM','MD','CB','++','+','--','-','TV'/
DATA COMMA,SPACE,NINE,ZERO/',',' ','9','0'/
TAB=CHAR(9)
IRCDRC=(INDEX(KEYWRD,'IRC')+INDEX(KEYWRD,'DRC') .NE.0)
ILOWA = ICHAR('a')
ILOWZ = ICHAR('z')
ICAPA = ICHAR('A')
ICAPZ = ICHAR('Z')
MAXTXT=0
NATOMS=0
NUMAT=0
ISERR=0
DO 10 I=1,MAXPAR
10 SIMBOL(I)= '---'
20 READ(IREAD,'(A)',END=130,ERR=230)LINE
IF(LINE.EQ.' ') GO TO 130
NATOMS=NATOMS+1
C
C SEE IF TEXT IS ASSOCIATED WITH THIS ELEMENT
C
I=INDEX(LINE,'(')
IF(I.NE.0)THEN
C
C YES, ELEMENT IS LABELLED.
C
K=INDEX(LINE,')')
TXTATM(NATOMS)=LINE(I:K)
MAXTXT=MAX(MAXTXT,K-I+1)
STRING=LINE(1:I-1)//LINE(K+1:)
LINE=STRING
ELSE
TXTATM(NATOMS)=' '
ENDIF
* CLEAN THE INPUT DATA
************************************************************************
DO 30 I=1,80
ILINE=ICHAR(LINE(I:I))
IF(ILINE.GE.ILOWA.AND.ILINE.LE.ILOWZ) THEN
LINE(I:I)=CHAR(ILINE+ICAPA-ILOWA)
ENDIF
30 CONTINUE
************************************************************************
ICOMMA=ICHAR(COMMA)
ITAB=ICHAR(TAB)
DO 40 I=1,80
KHAR=ICHAR(LINE(I:I))
IF(KHAR.EQ.ICOMMA.OR.KHAR.EQ.ITAB)LINE(I:I)=SPACE
40 CONTINUE
*
* INITIALIZE ISTART TO INTERPRET BLANKS AS ZERO'S
DO 50 I=1,10
50 ISTART(I)=80
*
* FIND INITIAL DIGIT OF ALL NUMBERS, CHECK FOR LEADING SPACES FOLLOWED
* BY A CHARACTER AND STORE IN ISTART
LEADSP=.TRUE.
NVALUE=0
DO 60 I=1,80
IF (LEADSP.AND.LINE(I:I).NE.SPACE) THEN
NVALUE=NVALUE+1
ISTART(NVALUE)=I
END IF
LEADSP=(LINE(I:I).EQ.SPACE)
60 CONTINUE
*
* ESTABLISH THE ELEMENT'S NAME AND ISOTOPE, CHECK FOR ERRORS OR E.O.DATA
*
WEIGHT=0.D0
STRING=LINE(ISTART(1):ISTART(2)-1)
IF( STRING(1:1) .GE. ZERO .AND. STRING(1:1) .LE. NINE) THEN
* ATOMIC NUMBER USED: NO ISOTOPE ALLOWED
LABEL=READA(STRING,1)
IF (LABEL.EQ.0) GO TO 120
IF (LABEL.LT.0.OR.LABEL.GT.107) THEN
WRITE(6,'('' ILLEGAL ATOMIC NUMBER'')')
GO TO 240
END IF
GO TO 80
END IF
* ATOMIC SYMBOL USED
REAL=ABS(READA(STRING,1))
IF (REAL.LT.1.D-15) THEN
* NO ISOTOPE
ELE=STRING(1:2)
ELSE
WEIGHT=REAL
IF( STRING(2:2) .GE. ZERO .AND. STRING(2:2) .LE. NINE) THEN
ELE=STRING(1:1)
ELSE
ELE=STRING(1:2)
END IF
END IF
* CHECK FOR ERROR IN ATOMIC SYMBOL
IF(ELE(1:1).EQ.'-'.AND.ELE(2:2).NE.'-')ELE(2:2)=' '
DO 70 I=1,107
IF(ELE.EQ.ELEMNT(I)) THEN
LABEL=I
GO TO 80
END IF
70 CONTINUE
IF(ELE(1:1).EQ.'X')THEN
LABEL=99
GOTO 80
ENDIF
WRITE(6,'('' UNRECOGNIZED ELEMENT NAME: ('',A,'')'')')ELE
GOTO 240
*
* ALL O.K.
*
80 IF (LABEL.NE.99) NUMAT=NUMAT+1
IF(WEIGHT.NE.0.D0)THEN
WRITE(6,'('' FOR ATOM'',I4,'' ISOTOPIC MASS:''
1 ,F15.5)')NATOMS, WEIGHT
ATMASS(NUMAT)=WEIGHT
ELSE
IF(LABEL .NE. 99) ATMASS(NUMAT)=AMS(LABEL)
ENDIF
IF(NATOMS.GT.NUMATM)THEN
WRITE(6,'(//10X,''**** MAX. NUMBER OF ATOMS ALLOWED:'',I4)')
1NUMATM
STOP
ENDIF
LABELS(NATOMS) =LABEL
GEO(1,NATOMS) =READA(LINE,ISTART(2))
GEO(2,NATOMS) =READA(LINE,ISTART(4))
GEO(3,NATOMS) =READA(LINE,ISTART(6))
IF(IRCDRC)THEN
TURN=LINE(ISTART(3):ISTART(3))
IF(TURN.EQ.'T')THEN
LOPT(1,NATOMS)=1
IF(NATOMS.EQ.1)WRITE(6,'(A)')' IN DRC MONITOR POTENTIAL ENER
1GY'//' TURNING POINTS'
ELSE
LOPT(1,NATOMS)=0
ENDIF
TURN=LINE(ISTART(5):ISTART(5))
IF(TURN.EQ.'T')THEN
LOPT(2,NATOMS)=1
ELSE
LOPT(2,NATOMS)=0
ENDIF
TURN=LINE(ISTART(7):ISTART(7))
IF(TURN.EQ.'T')THEN
LOPT(3,NATOMS)=1
ELSE
LOPT(3,NATOMS)=0
ENDIF
ELSE
LOPT(1,NATOMS) =READA(LINE,ISTART(3))
LOPT(2,NATOMS) =READA(LINE,ISTART(5))
LOPT(3,NATOMS) =READA(LINE,ISTART(7))
DO 90 I=3,7,2
IF(ICHAR(LINE(ISTART(I):ISTART(I))).GE.ICAPA.AND.
1ICHAR(LINE(ISTART(I):ISTART(I))).LE.ICAPZ)ISERR=1
90 CONTINUE
ENDIF
NA(NATOMS) =READA(LINE,ISTART(8))
NB(NATOMS) =READA(LINE,ISTART(9))
NC(NATOMS) =READA(LINE,ISTART(10))
C
C SPECIAL CASE OF USERS FORGETTING TO ADD DIHEDRAL DATA FOR ATOM 3
C
IF(NATOMS.EQ.3)THEN
IF(LOPT(3,3).EQ.2)THEN
NA(3)=1
NB(3)=2
GEO(3,3)=0.D0
LOPT(3,3)=0
ELSEIF(LOPT(3,3).EQ.1.AND.ABS(GEO(3,3)-2.D0).LT.1.D-4)THEN
NA(3)=2
NB(3)=1
GEO(3,3)=0.D0
LOPT(3,3)=0
ENDIF
ENDIF
IF(LOPT(1,NATOMS).GT.1.OR.LOPT(2,NATOMS).GT.1.OR.
1LOPT(3,NATOMS).GT.1)ISERR=1
IF(ISERR.EQ.1) THEN
C
C MUST BE GAUSSIAN GEOMETRY INPUT
C
DO 110 I=2,NATOMS
DO 110 K=1,3
J=GEO(K,I)+0.4D0
IF(ABS(GEO(K,I)-J).GT.1.D-5)THEN
C
C GEOMETRY CANNOT BE GAUSSIAN
C
WRITE(6,'(A)')' GEOMETRY IS FAULTY. GEOMETRY READ IN
1IS'
CONST=3.141592653598D0/180.D0
DO 100 L=1,NATOMS
GEO(2,L)=GEO(2,L)*CONST
100 GEO(3,L)=GEO(3,L)*CONST
CALL GEOUT(6)
STOP
ENDIF
110 CONTINUE
NATOMS=-1
RETURN
ENDIF
GOTO 20
*
* ALL DATA READ IN, CLEAN UP AND RETURN
*
120 NATOMS=NATOMS-1
130 NA(2)=1
LTXT=CHAR(MAXTXT)
IF(NATOMS.GT.3)THEN
INT=(NA(4).NE.0)
ELSE
IF(GEO(2,3).LT.10.AND.NATOMS.EQ.3)
1WRITE(6,'(//10X,'' WARNING: INTERNAL COORDINATES ARE ASSUMED -'',/
210X,'' FOR THREE-ATOM SYSTEMS '',//)')
INT=.TRUE.
ENDIF
IF(INT)GEO(2,2)=0
C
C READ IN VELOCITY VECTOR, IF PRESENT
C
IF(INDEX(KEYWRD,'VELO').GT.0)THEN
IF(INT)THEN
WRITE(6,'(A)')' COORDINATES MUST BE CARTESIAN WHEN VELOCITY'
1//' VECTOR IS USED.'
STOP
ENDIF
C# WRITE(6,'(/10X,A)')'INITIAL VELOCITY VECTOR FOR DRC'
DO 150 I=1,NATOMS
READ(5,'(A)') LINE
CALL NUCHAR(LINE,VALUE,NDMY)
IF(NDMY.NE.3)THEN
WRITE(6,'(/10X,A)')
1' THERE MUST BE EXACTLY THREE VELOCITY DATA PER LINE'
STOP
ENDIF
DO 140 J=1,3
140 REACT(J,I+2)=VALUE(J)
C# WRITE(6,'(2X,A2,2X,3F13.5)')ELEMNT(LABELS(I)),(VALUE(J),J=1,3)
150 CONTINUE
DO 160 I=1,3
DO 160 J=1,2
160 REACT(I,J)=GEO(I,J+1)-GEO(I,1)
C
C NOW TO ROTATE VELOCITY VECTOR TO SUIT INTERNAL COORDINATE DEFINITION
C
C
C ROTATE AROUND THE 1-2 X-AXIS TO AS TO ELIMINATE REACT(3,2)
C (PUT ATOM 2 IN X-Y PLANE)
SA=REACT(3,1)/SQRT(REACT(2,1)**2+REACT(3,1)**2+1.D-20)
CA=SIGN(SQRT(1.D0-SA**2),REACT(2,1))
C# LABELS(NATOMS+1)=1
C# LABELS(NATOMS+2)=1
C# WRITE(6,*)' FIRST ROTATION, ABOUT 1-2 X-AXIS'
DO 170 I=1,NATOMS+2
TEMP1= REACT(2,I)*CA+REACT(3,I)*SA
TEMP2=-REACT(2,I)*SA+REACT(3,I)*CA
REACT(2,I)=TEMP1
REACT(3,I)=TEMP2
C# WRITE(6,'(2X,A2,2X,3F13.5)')ELEMNT(LABELS(I)),(REACT(J,I),J=1,3)
170 CONTINUE
C ROTATE AROUND THE 1-2 Z-AXIS TO AS TO ELIMINATE REACT(2,2)
C (PUT ATOM 2 ON X AXIS)
CA=REACT(1,1)/SQRT(REACT(2,1)**2+REACT(1,1)**2+1.D-20)
SA=SIGN(SQRT(1.D0-CA**2),REACT(2,1))
C# WRITE(6,*)' SECOND ROTATION, ABOUT 1-2 Z-AXIS'
DO 180 I=1,NATOMS+2
TEMP1= REACT(1,I)*CA+REACT(2,I)*SA
TEMP2=-REACT(1,I)*SA+REACT(2,I)*CA
REACT(1,I)=TEMP1
REACT(2,I)=TEMP2
C# WRITE(6,'(2X,A2,2X,3F13.5)')ELEMNT(LABELS(I)),(REACT(J,I),J=1,3)
180 CONTINUE
C ROTATE AROUND THE 2-3 X-AXIS TO AS TO ELIMINATE REACT(3,3)
C (PUT ATOM 3 ON X-Y PLANE)
SA=REACT(3,2)/SQRT(REACT(2,2)**2+REACT(3,2)**2+1.D-20)
CA=SIGN(SQRT(1.D0-SA**2),REACT(2,2))
C# WRITE(6,*)' THIRD ROTATION, ABOUT 2-3 X-AXIS'
DO 190 I=1,NATOMS+2
TEMP1= REACT(2,I)*CA+REACT(3,I)*SA
TEMP2=-REACT(2,I)*SA+REACT(3,I)*CA
REACT(2,I)=TEMP1
REACT(3,I)=TEMP2
C# WRITE(6,'(2X,A2,2X,3F13.5)')ELEMNT(LABELS(I)),(REACT(J,I),J=1,3)
190 CONTINUE
C
C STRIP OFF FIRST TWO COORDINATES; THESE WERE THE COORDINATE AXIS
C DEFINITIONS
C
DO 200 I=1,NATOMS
DO 200 J=1,3
200 REACT(J,I)=REACT(J,I+2)
ENDIF
IF( .NOT. INT ) THEN
DO 210 I=1,NATOMS
DO 210 J=1,3
210 XYZ(J,I)=GEO(J,I)
DEGREE=90.D0/ASIN(1.D0)
CALL XYZINT(XYZ,NATOMS,NA,NB,NC,DEGREE,GEO)
IF(INDEX(KEYWRD,' XYZ').EQ.0)THEN
C
C UNCONDITIONALLY SET FLAGS FOR INTERNAL COORDINATES
C
DO 220 I=1,3
DO 220 J=I,3
220 LOPT(J,I)=0
ENDIF
IF(ABS(GEO(2,3)-180.D0).LT.1.D-4.OR.ABS(GEO(2,3)).LT.1.D-4)
1THEN
WRITE(6,'(A)')' DUE TO PROGRAM BUG, THE FIRST THREE ATOMS MU
1ST NOT LIE IN A STRAIGHT LINE.'
STOP
ENDIF
ELSEIF (.NOT.IRCDRC) THEN
LOPT(2,2)=0
IF(LOPT(1,1)+LOPT(2,1)+LOPT(3,1)+LOPT(3,2)+
1 LOPT(3,3) .GT. 0)THEN
LOPT(1,1)=0
LOPT(2,1)=0
LOPT(3,1)=0
LOPT(3,2)=0
LOPT(3,3)=0
WRITE(6,'(//10X,'' AN UNOPTIMIZABLE GEOMETRIC PARAMETER HAS'
1',/10X,'' BEEN MARKED FOR OPTIMIZATION. THIS IS A NON-FATAL ''
2,''ERROR'')')
ENDIF
ENDIF
IF(NA(3).EQ.0) THEN
NB(3)=1
NA(3)=2
ENDIF
RETURN
* ERROR CONDITIONS
230 IF(IREAD.EQ.5) THEN
WRITE(6,'( '' ERROR DURING READ AT ATOM NUMBER '', I3 )')NATOMS
ELSE
NATOMS=0
RETURN
ENDIF
240 J=NATOMS-1
WRITE(6,'('' DATA CURRENTLY READ IN ARE '')')
DO 250 K=1,J
250 WRITE(6,260)LABELS(K),(GEO(JJ,K),LOPT(JJ,K),JJ=1,3),
1NA(K),NB(K),NC(K)
260 FORMAT(I4,2X,3(F10.5,2X,I2,2X),3(I2,1X))
STOP
END
|
