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
|
SUBROUTINE XYZINT(XYZ,NUMAT,NA,NB,NC,DEGREE,GEO)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION XYZ(3,*), NA(*), NB(*), NC(*), GEO(3,*)
***********************************************************************
*
* XYZINT WORKS OUT THE INTERNAL COORDINATES OF A MOLECULE.
* THE "RULES" FOR THE CONNECTIVITY ARE AS FOLLOWS:
* ATOM I IS DEFINED AS BEING AT A DISTANCE FROM THE NEAREST
* ATOM J, ATOM J ALREADY HAVING BEEN DEFINED.
* ATOM I MAKES AN ANGLE WITH ATOM J AND THE ATOM K, WHICH HAS
* ALREADY BEEN DEFINED, AND IS THE NEAREST ATOM TO J
* ATOM I MAKES A DIHEDRAL ANGLE WITH ATOMS J, K, AND L. L HAVING
* BEEN DEFINED AND IS THE NEAREST ATOM TO K, AND J, K AND L
* HAVE A CONTAINED ANGLE IN THE RANGE 15 TO 165 DEGREES,
* IF POSSIBLE.
*
* IF(NA(2).EQ.-1 OR -2 THEN THE ORIGINAL CONNECTIVITY IS USED.
*
* NOTE THAT GEO AND XYZ MUST NOT BE THE SAME IN THE CALL.
*
* ON INPUT XYZ = CARTESIAN ARRAY OF NUMAT ATOMS
* DEGREE = 1 IF ANGLES ARE TO BE IN RADIANS
* DEGREE = 57.29578 IF ANGLES ARE TO BE IN DEGREES
*
***********************************************************************
COMMON /GEOOK/ IGEOOK
COMMON /NUMCAL/ NUMCAL
DATA ICALCN/0/
IGEOOK=99
IF(.NOT.(ICALCN.NE.NUMCAL).AND.NA(2).EQ.-1 .OR. NA(2).EQ.-2)THEN
NA(2)=1
DO 10 I=2,NUMAT
J=NA(I)
IF(I.GT.3)CALL DIHED(XYZ,I,J,NB(I),NC(I),GEO(3,I))
IF(I.GT.2)CALL BANGLE(XYZ,I,J,NB(I),GEO(2,I))
GEO(1,I)=SQRT((XYZ(1,I)-XYZ(1,J))**2+
1 (XYZ(2,I)-XYZ(2,J))**2+
2 (XYZ(3,I)-XYZ(3,J))**2)
10 CONTINUE
ELSE
IF(NA(2).EQ.-1)ICALCN=NUMCAL
DO 30 I=1,NUMAT
NA(I)=2
NB(I)=3
NC(I)=4
IM1=I-1
IF(IM1.EQ.0)GOTO 30
SUM=1.D30
DO 20 J=1,IM1
R=(XYZ(1,I)-XYZ(1,J))**2+
1 (XYZ(2,I)-XYZ(2,J))**2+
2 (XYZ(3,I)-XYZ(3,J))**2
IF(R.LT.SUM.AND.NA(J).NE.J.AND.NB(J).NE.J) THEN
SUM=R
K=J
ENDIF
20 CONTINUE
C
C ATOM I IS NEAREST TO ATOM K
C
NA(I)=K
IF(I.GT.2)NB(I)=NA(K)
IF(I.GT.3)NC(I)=NB(K)
C
C FIND ANY ATOM TO RELATE TO NA(I)
C
30 CONTINUE
ENDIF
NA(1)=0
NB(1)=0
NC(1)=0
NB(2)=0
NC(2)=0
NC(3)=0
CALL XYZGEO(XYZ,NUMAT,NA,NB,NC,DEGREE,GEO)
RETURN
END
SUBROUTINE XYZGEO(XYZ,NUMAT,NA,NB,NC,DEGREE,GEO)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION XYZ(3,*), NA(*), NB(*), NC(*), GEO(3,*)
***********************************************************************
*
* XYZGEO CONVERTS COORDINATES FROM CARTESIAN TO INTERNAL.
*
* ON INPUT XYZ = ARRAY OF CARTESIAN COORDINATES
* NUMAT= NUMBER OF ATOMS
* NA = NUMBERS OF ATOM TO WHICH ATOMS ARE RELATED
* BY DISTANCE
* NB = NUMBERS OF ATOM TO WHICH ATOMS ARE RELATED
* BY ANGLE
* NC = NUMBERS OF ATOM TO WHICH ATOMS ARE RELATED
* BY DIHEDRAL
*
* ON OUTPUT GEO = INTERNAL COORDINATES IN ANGSTROMS, RADIANS,
* AND RADIANS
*
***********************************************************************
DO 30 I=2,NUMAT
J=NA(I)
K=NB(I)
L=NC(I)
IF(I.LT.3) GOTO 30
II=I
CALL BANGLE(XYZ,II,J,K,GEO(2,I))
GEO(2,I)=GEO(2,I)*DEGREE
IF(I.LT.4) GOTO 30
C
C MAKE SURE DIHEDRAL IS MEANINGLFUL
C
CALL BANGLE(XYZ,J,K,L,ANGL)
TOL=0.2617994D0
IF(ANGL.GT.3.1415926D0-TOL.OR.ANGL.LT.TOL)THEN
C
C ANGLE IS UNSATISFACTORY, LET'S SEARCH FOR ANOTHER ATOM FOR
C DEFINING THE DIHEDRAL.
10 SUM=100.D0
DO 20 I1=1,II-1
R=(XYZ(1,I1)-XYZ(1,K))**2+
1 (XYZ(2,I1)-XYZ(2,K))**2+
2 (XYZ(3,I1)-XYZ(3,K))**2
IF(R.LT.SUM.AND.I1.NE.J.AND.I1.NE.K) THEN
CALL BANGLE(XYZ,J,K,I1,ANGL)
IF(ANGL.LT.3.1415926D0-TOL.AND.ANGL.GT.TOL)THEN
SUM=R
L=I1
NC(II)=L
ENDIF
ENDIF
20 CONTINUE
IF(SUM.GT.99.D0.AND.TOL.GT.0.1D0)THEN
C
C ANYTHING WITHIN 5 DEGREES?
C
TOL=0.087266D0
GOTO 10
ENDIF
ENDIF
CALL DIHED(XYZ,II,J,K,L,GEO(3,I))
GEO(3,I)=GEO(3,I)*DEGREE
30 GEO(1,I)= SQRT((XYZ(1,I)-XYZ(1,J))**2+
1 (XYZ(2,I)-XYZ(2,J))**2+
2 (XYZ(3,I)-XYZ(3,J))**2)
GEO(1,1)=0.D0
GEO(2,1)=0.D0
GEO(3,1)=0.D0
GEO(2,2)=0.D0
GEO(3,2)=0.D0
GEO(3,3)=0.D0
RETURN
END
SUBROUTINE BANGLE(XYZ,I,J,K,ANGLE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION XYZ(3,*)
*********************************************************************
*
* BANGLE CALCULATES THE ANGLE BETWEEN ATOMS I,J, AND K. THE
* CARTESIAN COORDINATES ARE IN XYZ.
*
*********************************************************************
D2IJ = (XYZ(1,I)-XYZ(1,J))**2+
1 (XYZ(2,I)-XYZ(2,J))**2+
2 (XYZ(3,I)-XYZ(3,J))**2
D2JK = (XYZ(1,J)-XYZ(1,K))**2+
1 (XYZ(2,J)-XYZ(2,K))**2+
2 (XYZ(3,J)-XYZ(3,K))**2
D2IK = (XYZ(1,I)-XYZ(1,K))**2+
1 (XYZ(2,I)-XYZ(2,K))**2+
2 (XYZ(3,I)-XYZ(3,K))**2
XY = SQRT(D2IJ*D2JK)
TEMP = 0.5D0 * (D2IJ+D2JK-D2IK) / XY
IF (TEMP .GT. 1.0D0) TEMP=1.0D0
IF (TEMP .LT. -1.0D0) TEMP=-1.0D0
ANGLE = ACOS( TEMP )
RETURN
END
SUBROUTINE DIHED(XYZ,I,J,K,L,ANGLE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION XYZ(3,*)
*********************************************************************
*
* DIHED CALCULATES THE DIHEDRAL ANGLE BETWEEN ATOMS I, J, K,
* AND L. THE CARTESIAN COORDINATES OF THESE ATOMS
* ARE IN ARRAY XYZ.
*
* DIHED IS A MODIFIED VERSION OF A SUBROUTINE OF THE SAME NAME
* WHICH WAS WRITTEN BY DR. W. THEIL IN 1973.
*
*********************************************************************
XI1=XYZ(1,I)-XYZ(1,K)
XJ1=XYZ(1,J)-XYZ(1,K)
XL1=XYZ(1,L)-XYZ(1,K)
YI1=XYZ(2,I)-XYZ(2,K)
YJ1=XYZ(2,J)-XYZ(2,K)
YL1=XYZ(2,L)-XYZ(2,K)
ZI1=XYZ(3,I)-XYZ(3,K)
ZJ1=XYZ(3,J)-XYZ(3,K)
ZL1=XYZ(3,L)-XYZ(3,K)
C ROTATE AROUND Z AXIS TO PUT KJ ALONG Y AXIS
DIST= SQRT(XJ1**2+YJ1**2+ZJ1**2)
COSA=ZJ1/DIST
IF(COSA.GT.1.0D0) COSA=1.0D0
IF(COSA.LT.-1.0D0) COSA=-1.0D0
DDD=1.0D0-COSA**2
IF(DDD.LE.0.0) GO TO 10
YXDIST=DIST* SQRT(DDD)
IF(YXDIST.GT.1.0D-6) GO TO 20
10 CONTINUE
XI2=XI1
XL2=XL1
YI2=YI1
YL2=YL1
COSTH=COSA
SINTH=0.D0
GO TO 30
20 COSPH=YJ1/YXDIST
SINPH=XJ1/YXDIST
XI2=XI1*COSPH-YI1*SINPH
XL2=XL1*COSPH-YL1*SINPH
YI2=XI1*SINPH+YI1*COSPH
YJ2=XJ1*SINPH+YJ1*COSPH
YL2=XL1*SINPH+YL1*COSPH
C ROTATE KJ AROUND THE X AXIS SO KJ LIES ALONG THE Z AXIS
COSTH=COSA
SINTH=YJ2/DIST
30 CONTINUE
YI3=YI2*COSTH-ZI1*SINTH
YL3=YL2*COSTH-ZL1*SINTH
CALL DANG(XL2,YL3,XI2,YI3,ANGLE)
IF (ANGLE .LT. 0.) ANGLE=4.0D0* ASIN(1.0D00)+ANGLE
IF (ANGLE .GE. 6.2831853D0 ) ANGLE=0.D0
RETURN
END
SUBROUTINE DANG(A1,A2,B1,B2,RCOS)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
**********************************************************************
*
* DANG DETERMINES THE ANGLE BETWEEN THE POINTS (A1,A2), (0,0),
* AND (B1,B2). THE RESULT IS PUT IN RCOS.
*
**********************************************************************
ZERO=1.0D-6
IF( ABS(A1).LT.ZERO.AND. ABS(A2).LT.ZERO) GO TO 10
IF( ABS(B1).LT.ZERO.AND. ABS(B2).LT.ZERO) GO TO 10
ANORM=1.0D0/ SQRT(A1**2+A2**2)
BNORM=1.0D0/ SQRT(B1**2+B2**2)
A1=A1*ANORM
A2=A2*ANORM
B1=B1*BNORM
B2=B2*BNORM
SINTH=(A1*B2)-(A2*B1)
COSTH=A1*B1+A2*B2
IF(COSTH.GT.1.0D0) COSTH=1.0D0
IF(COSTH.LT.-1.0D0) COSTH=-1.0D0
RCOS= ACOS(COSTH)
IF( ABS(RCOS).LT.4.0D-4) GO TO 10
IF(SINTH.GT.0.D0) RCOS=4.0D0* ASIN(1.0D00)-RCOS
RCOS=-RCOS
RETURN
10 RCOS=0.0D0
RETURN
END
|
