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
|
SUBROUTINE GMETRY(GEO,COORD)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'SIZES'
COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
1NA(NUMATM),NB(NUMATM),NC(NUMATM)
2 /EULER / TVEC(3,3), ID
COMMON /REACTN/ STEP, GEOA(3,NUMATM), GEOVEC(3,NUMATM),COLCST
COMMON /GEOOK/ IGEOOK
COMMON /NUMCAL/ NUMCAL
DIMENSION GEO(3,*),COORD(3,*)
CHARACTER *15 NDIMEN(4)
LOGICAL GEOOK
SAVE ICALCN, NDIMEN, GEOOK
DATA ICALCN/0/
DATA NDIMEN/' MOLECULE ',' POLYMER ',
1'LAYER STRUCTURE',' SOLID '/
C***********************************************************************
C
C GMETRY COMPUTES COORDINATES FROM BOND-ANGLES AND LENGTHS.
C *** IT IS ADAPTED FROM THE PROGRAM WRITTEN BY M.J.S. DEWAR.
C
C THREE SEPARATE OPTIONS EXIST WITHIN GMETRY. THESE ARE:
C (A) IF NA(1) IS EQUAL TO 99 (IMPOSSIBLE UNDER NORMAL CIRCUMSTANCES)
C THEN GEO IS ASSUMED TO BE IN CARTESIAN RATHER THAN INTERNAL
C COORDINATES, AND COORD IS THEN SET EQUAL TO GEO.
C (B) IF STEP IS NON-ZERO (THIS IS THE CASE WHEN "SADDLE" IS USED)
C THEN GEO IS FIRST MODIFIED BY SHIFTING THE INTERNAL COORDINATES
C ALONG A RADIUS FROM GEOA TO PLACE GEO AT ADISTANCESTEPFROMGEOA.
C (C) NORMAL CONVERSION FROM INTERNAL TO CARTESIAN COORDINATESISDONE.
C
C ON INPUT:
C GEO = ARRAY OF INTERNAL COORDINATES.
C NATOMS = NUMBER OF ATOMS, INCLUDING DUMMIES.
C NA = ARRAY OF ATOM LABELS FOR BOND LENGTHS.
C
C ON OUTPUT:
C COORD = ARRAY OF CARTESIAN COORDINATES
C
C***********************************************************************
C OPTION (B)
GEOOK=(IGEOOK.EQ.99)
IF(ABS(STEP) .GT. 1.D-4) THEN
SUM=0.D0
DO 10 J=1,3
C$DOIT VBEST
DO 10 I=1,NATOMS
GEOVEC(J,I)=GEO(J,I)-GEOA(J,I)
10 SUM=SUM+GEOVEC(J,I)**2
SUM=SQRT(SUM)
ERROR=(SUM-STEP)/SUM
ELSE
ERROR=0.D0
ENDIF
DO 20 J=1,3
C$DOIT VBEST
DO 20 I=1,NATOMS
20 GEO(J,I)=GEO(J,I)-ERROR*GEOVEC(J,I)
C OPTION (A)
IF(NA(1).EQ.99) THEN
DO 30 I=1,3
C$DOIT VBEST
DO 30 J=1,NATOMS
30 COORD(I,J)=GEO(I,J)
GOTO 100
ENDIF
C OPTION (C)
COORD(1,1)=0.0D00
COORD(2,1)=0.0D00
COORD(3,1)=0.0D00
COORD(1,2)=GEO(1,2)
COORD(2,2)=0.0D00
COORD(3,2)=0.0D00
IF(NATOMS.EQ.2) GOTO 100
CCOS=COS(GEO(2,3))
IF(NA(3).EQ.1)THEN
COORD(1,3)=COORD(1,1)+GEO(1,3)*CCOS
ELSE
COORD(1,3)=COORD(1,2)-GEO(1,3)*CCOS
ENDIF
COORD(2,3)=GEO(1,3)*SIN(GEO(2,3))
COORD(3,3)=0.0D00
DO 90 I=4,NATOMS
COSA=COS(GEO(2,I))
MB=NB(I)
MC=NA(I)
XB=COORD(1,MB)-COORD(1,MC)
YB=COORD(2,MB)-COORD(2,MC)
ZB=COORD(3,MB)-COORD(3,MC)
RBC=XB*XB+YB*YB+ZB*ZB
IF(RBC.LT.1.D-16)THEN
C
C TWO ATOMS ARE COINCIDENT. A FATAL ERROR.
C
WRITE(6,'(A,I4,A,I4,A)')' ATOMS',MB,' AND',MC,' ARE COINCIDE
1NT'
WRITE(6,'(A)')' THIS IS A FATAL ERROR, RUN STOPPED IN GMETRY
1'
STOP
ELSE
RBC=1.0D00/SQRT(RBC)
ENDIF
MA=NC(I)
XA=COORD(1,MA)-COORD(1,MC)
YA=COORD(2,MA)-COORD(2,MC)
ZA=COORD(3,MA)-COORD(3,MC)
C
C ROTATE ABOUT THE Z-AXIS TO MAKE YB=0, AND XB POSITIVE. IF XYB IS
C TOO SMALL, FIRST ROTATE THE Y-AXIS BY 90 DEGREES.
C
XYB=SQRT(XB*XB+YB*YB)
K=-1
IF (XYB.GT.0.1D00) GO TO 40
XPA=ZA
ZA=-XA
XA=XPA
XPB=ZB
ZB=-XB
XB=XPB
XYB=SQRT(XB*XB+YB*YB)
K=+1
C
C ROTATE ABOUT THE Y-AXIS TO MAKE ZB VANISH
C
40 COSTH=XB/XYB
SINTH=YB/XYB
XPA=XA*COSTH+YA*SINTH
YPA=YA*COSTH-XA*SINTH
SINPH=ZB*RBC
COSPH=SQRT(ABS(1.D00-SINPH*SINPH))
ZQA=ZA*COSPH-XPA*SINPH
C
C ROTATE ABOUT THE X-AXIS TO MAKE ZA=0, AND YA POSITIVE.
C
YZA=SQRT(YPA**2+ZQA**2)
IF(YZA.LT.1.D-4)GOTO 60
IF(YZA.LT.2.D-2 .AND. .NOT.GEOOK)THEN
WRITE(6,'(//20X,'' CALCULATION ABANDONED AT THIS POINT'')')
WRITE(6,'(//10X,'' THREE ATOMS BEING USED TO DEFINE THE'',/
110X,'' COORDINATES OF A FOURTH ATOM, WHOSE BOND-ANGLE IS'')')
WRITE(6,'(10X,'' NOT ZERO OR 180 DEGREEES, ARE '',
1''IN AN ALMOST STRAIGHT'')')
WRITE(6,'(10X,'' LINE. THERE IS A HIGH PROBABILITY THAT THE
1'',/10X,'' COORDINATES OF THE ATOM WILL BE INCORRECT.'')')
WRITE(6,'(//20X,''THE FAULTY ATOM IS ATOM NUMBER'',I4)')I
CALL GEOUT(1)
WRITE(6,'(//20X,''CARTESIAN COORDINATES UP TO FAULTY ATOM'')
1')
WRITE(6,'(//5X,''I'',12X,''X'',12X,''Y'',12X,''Z'')')
DO 50 J=1,I
50 WRITE(6,'(I6,F16.5,2F13.5)')J,(COORD(K,J),K=1,3)
WRITE(6,'(//6X,'' ATOMS'',I3,'','',I3,'', AND'',I3,
1'' ARE WITHIN'',F7.4,'' ANGSTROMS OF A STRAIGHT LINE'')')
2MC,MB,MA,YZA
STOP
ENDIF
COSKH=YPA/YZA
SINKH=ZQA/YZA
GOTO 70
60 CONTINUE
C
C ANGLE TOO SMALL TO BE IMPORTANT
C
COSKH=1.D0
SINKH=0.D0
70 CONTINUE
C
C COORDINATES :- A=(???,YZA,0), B=(RBC,0,0), C=(0,0,0)
C NONE ARE NEGATIVE.
C THE COORDINATES OF I ARE EVALUATED IN THE NEW FRAME.
C
SINA=SIN(GEO(2,I))
SIND=-SIN(GEO(3,I))
COSD=COS(GEO(3,I))
XD=GEO(1,I)*COSA
YD=GEO(1,I)*SINA*COSD
ZD=GEO(1,I)*SINA*SIND
C
C TRANSFORM THE COORDINATES BACK TO THE ORIGINAL SYSTEM.
C
YPD=YD*COSKH-ZD*SINKH
ZPD=ZD*COSKH+YD*SINKH
XPD=XD*COSPH-ZPD*SINPH
ZQD=ZPD*COSPH+XD*SINPH
XQD=XPD*COSTH-YPD*SINTH
YQD=YPD*COSTH+XPD*SINTH
IF (K.LT.1) GO TO 80
XRD=-ZQD
ZQD=XQD
XQD=XRD
80 COORD(1,I)=XQD+COORD(1,MC)
COORD(2,I)=YQD+COORD(2,MC)
COORD(3,I)=ZQD+COORD(3,MC)
90 CONTINUE
C
C *** NOW REMOVE THE TRANSLATION VECTORS, IF ANY, FROM THE ARRAY COOR
C
100 CONTINUE
K=NATOMS
110 IF(LABELS(K).NE.107) GOTO 120
K=K-1
GOTO 110
120 K=K+1
IF(K.GT.NATOMS) GOTO 170
C
C SYSTEM IS A SOLID, OF DIMENSION NATOMS+1-K
C
L=0
DO 130 I=K,NATOMS
L=L+1
MC=NA(I)
TVEC(1,L)=COORD(1,I)-COORD(1,MC)
TVEC(2,L)=COORD(2,I)-COORD(2,MC)
TVEC(3,L)=COORD(3,I)-COORD(3,MC)
130 CONTINUE
ID=L
IF (ICALCN.NE.NUMCAL) THEN
ICALCN=NUMCAL
WRITE(6,140)NDIMEN(ID+1)
140 FORMAT(/10X,' THE SYSTEM IS A ',A15,/)
IF(ID.EQ.0) GOTO 170
WRITE(6,150)
WRITE(6,160)(I,(TVEC(J,I),J=1,3),I=1,ID)
150 FORMAT(/,' UNIT CELL TRANSLATION VECTORS',/
1/,' X Y Z')
160 FORMAT(' T',I1,' = ',F11.7,' ',F11.7,' ',F11.7)
ENDIF
170 CONTINUE
J=0
DO 190 I=1,NATOMS
IF (LABELS(I).EQ.99.OR.LABELS(I).EQ.107) GO TO 190
J=J+1
C$DOIT ASIS
DO 180 K=1,3
180 COORD(K,J)=COORD(K,I)
190 CONTINUE
RETURN
END
|
