1 |
#define EXTERN extern |
2 |
#include "pcwin.h" |
3 |
#include "pcmod.h" |
4 |
#include "field.h" |
5 |
#include "atom_k.h" |
6 |
#include "draw.h" |
7 |
|
8 |
void hadd(void); |
9 |
void hdel(int); |
10 |
void type(void); |
11 |
void nhadd(void); |
12 |
void hcoord(void); |
13 |
static void revec(void); |
14 |
static void xyplan(float *, float *, float *); |
15 |
static void xaxis(float *,float *,float *); |
16 |
static void getvec(void); |
17 |
static void reseq(void); |
18 |
|
19 |
EXTERN struct t_minim_control { |
20 |
int type, method, field, added_const; |
21 |
char added_path[256],added_name[256]; |
22 |
} minim_control; |
23 |
|
24 |
void hadd(void) |
25 |
{ |
26 |
type(); |
27 |
if (field.type == MM3) |
28 |
{ |
29 |
nhadd(); |
30 |
hcoord(); |
31 |
hdel(1); |
32 |
} else if (field.type == MMFF94) |
33 |
{ |
34 |
nhadd(); |
35 |
hcoord(); |
36 |
hdel(1); |
37 |
} else |
38 |
{ |
39 |
set_atomtypes(MMX); |
40 |
nhadd(); |
41 |
hcoord(); |
42 |
} |
43 |
type(); |
44 |
} |
45 |
// ======================== |
46 |
void nhadd() |
47 |
{ |
48 |
int i, i1, ik, im, it, itads, j, jk, nhadds, newatom; |
49 |
int jji, bo,ncount; |
50 |
static int iadd[MAXVDW]={ 0,0,0,0,0,3,1,6,0,1, |
51 |
0,0,0,0,3,6,3,0,0,0, |
52 |
0,0,0,0,6,0,0,0,0,0, |
53 |
0,0,0,3,0,0,6,1,1,0, |
54 |
0,4,0,0,0,6,0,6,0,0, |
55 |
3,3,3,0,6,0,0,0,0,0, |
56 |
0,0,0,0,0,1,6,0,0,0, |
57 |
0,0,0,0,0}; |
58 |
|
59 |
/* this routine calculates the #'s of h to be added |
60 |
* to carbon and packs it symbolically into iat(i,j) table |
61 |
* as well as into the connection table via the appropriate |
62 |
* manipulation instructions. hydrogen coordinates |
63 |
* are not added by this routine. lp are also added to oxygen. */ |
64 |
/* 0 to disregard : 1 and 4 for lp's only |
65 |
* 2 and 5 for h's only and 3 and 6 for h and lp |
66 |
* add 3 for trigonal atoms */ |
67 |
|
68 |
ncount = natom; |
69 |
for( i = 1 ; i <= ncount; i++ ) |
70 |
{ |
71 |
if( atom[i].mmx_type > 0 && (atom[i].mmx_type < 300) && atom[i].atomnum != 1 ) |
72 |
{ |
73 |
jji = 0; |
74 |
bo = 0; |
75 |
for (j=0; j < MAXIAT; j++) |
76 |
{ |
77 |
if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9) |
78 |
{ |
79 |
jji ++; |
80 |
bo += atom[i].bo[j]; |
81 |
} |
82 |
} |
83 |
nhadds = atom_def.ligands[atom[i].mmx_type] - jji; |
84 |
|
85 |
if (atom[i].mmx_type >= 11 && atom[i].mmx_type <= 14) |
86 |
nhadds = 0; |
87 |
if (atom[i].mmx_type == 41 || atom[i].mmx_type == 46) |
88 |
nhadds = atom_def.ligands[atom[i].mmx_type] - bo; |
89 |
if (atom[i].mmx_type == 49 || atom[i].mmx_type == 50 || atom[i].mmx_type == 51) |
90 |
nhadds = atom_def.ligands[atom[i].mmx_type] - bo; |
91 |
|
92 |
/* add atom and bond it to atom i */ |
93 |
if( nhadds > 0 ) |
94 |
{ |
95 |
for( im = 1; im <= nhadds; im++ ) |
96 |
{ |
97 |
newatom = make_atom(5,0.F,0.F,0.F,"H"); |
98 |
set_atomdata(newatom,5,5,5); |
99 |
if (newatom == -1) |
100 |
{ |
101 |
message_alert("Error in nhadd newatom = -1","Error"); |
102 |
hdel(0); |
103 |
return; |
104 |
} |
105 |
make_bond(i,newatom,1); |
106 |
} |
107 |
} |
108 |
} |
109 |
} |
110 |
|
111 |
/* *** check for lp's *** */ |
112 |
if (field.type == MMX) |
113 |
{ |
114 |
ncount = natom; |
115 |
for( i = 1 ; i <= ncount; i++ ) |
116 |
{ |
117 |
it = iadd[atom[i].mmx_type-1]; |
118 |
if( it >= 1 ) |
119 |
{ |
120 |
itads = atom[i].mmx_type ; |
121 |
if (it == 1 || it == 3 || it == 4 || it == 6 ) |
122 |
{ |
123 |
nhadds = 4; |
124 |
for( i1 = 0; i1 < MAXIAT; i1++ ) |
125 |
{ |
126 |
if( atom[i].bo[i1] != 0 && atom[i].bo[i1] != 9 ) |
127 |
nhadds = nhadds - atom[i].bo[i1]; |
128 |
} |
129 |
if( itads == 6 || itads == 15 || itads == 34 || itads == 42 || itads == 53 || itads == 48 ) |
130 |
{ |
131 |
for( ik = 0; ik < MAXIAT; ik++ ) |
132 |
{ |
133 |
if( atom[i].iat[ik] != 0 && atom[i].bo[ik] != 9 ) |
134 |
{ |
135 |
jk = atom[atom[i].iat[ik]].mmx_type; |
136 |
if( jk == 2 || jk == 3 || jk == 4 || jk == 9 |
137 |
|| jk == 10 || jk == 30 || jk == 37 || jk == 40 ) |
138 |
{ |
139 |
nhadds = nhadds - 1; |
140 |
break; |
141 |
} |
142 |
} |
143 |
} |
144 |
} |
145 |
|
146 |
/* correction for sulfoxide */ |
147 |
if( itads == 17 ) |
148 |
nhadds = nhadds + 1; |
149 |
|
150 |
if (itads == 66) |
151 |
nhadds = 2; |
152 |
|
153 |
if( nhadds > 0 ) |
154 |
{ |
155 |
for( i1 = 1; i1 <= nhadds; i1++ ) |
156 |
{ |
157 |
newatom = make_atom(20, 0.F,0.F,0.F,""); |
158 |
set_atomdata(newatom,20,0,0); |
159 |
make_bond(i,newatom,1); |
160 |
} |
161 |
} |
162 |
} |
163 |
} |
164 |
} |
165 |
} |
166 |
return; |
167 |
} |
168 |
// ================================== |
169 |
void hcoord() |
170 |
{ |
171 |
int i, i1, i3, i4, i5, i6, i7, i8, |
172 |
i9, iadj[10], ii, ii1, ii2, ii3, iii, |
173 |
ij, ijk[3], it, j, j1, j2, j3, j4, j5, j8, |
174 |
j9, k1, k2, l1, na, nhadds, ni; |
175 |
float adj, dis, xpnt, ypnt, zpnt; |
176 |
|
177 |
getvec(); |
178 |
for( i = 1; i <= natom; i++ ) |
179 |
{ |
180 |
it = atom[i].mmx_type; |
181 |
if( it < 300) |
182 |
{ |
183 |
for( j = 0; j < 10; j++ ) |
184 |
iadj[j] = 0; |
185 |
|
186 |
for( j = 0; j < 3; j++ ) |
187 |
ijk[j] = 0; |
188 |
|
189 |
nhadds = 0; |
190 |
l1 = 0; |
191 |
for( j = 0; j < MAXIAT; j++ ) |
192 |
{ |
193 |
if( atom[i].iat[j] != 0 && atom[i].bo[j] != 9 ) |
194 |
{ |
195 |
if( atom[atom[i].iat[j]].mmx_type == 5 || atom[atom[i].iat[j]].mmx_type == 20 ) |
196 |
nhadds++ ; |
197 |
else |
198 |
{ |
199 |
iadj[l1] = atom[i].iat[j]; |
200 |
l1++; |
201 |
} |
202 |
} |
203 |
} |
204 |
if( nhadds != 0 ) |
205 |
{ |
206 |
/* **** calculate the coordinates of the center we |
207 |
* are going to move to the x-axis */ |
208 |
xpnt = atom[i].x; |
209 |
ypnt = atom[i].y; |
210 |
zpnt = atom[i].z; |
211 |
for( ii = 1; ii <= (natom + 4); ii++ ) |
212 |
{ |
213 |
atom[ii].x = atom[ii].x - xpnt; |
214 |
atom[ii].y = atom[ii].y - ypnt; |
215 |
atom[ii].z = atom[ii].z - zpnt; |
216 |
} |
217 |
xpnt = 0.F; |
218 |
ypnt = 0.F; |
219 |
zpnt = 0.F; |
220 |
adj = 0.F; |
221 |
for( i1 = 0; i1 < 10; i1++ ) |
222 |
{ |
223 |
if( iadj[i1] != 0 ) |
224 |
{ |
225 |
na = iadj[i1]; |
226 |
adj += 1.0001F ; |
227 |
xpnt += atom[na].x; |
228 |
ypnt += atom[na].y; |
229 |
zpnt += atom[na].z; |
230 |
} |
231 |
} |
232 |
if( adj > 0.01F ) |
233 |
{ |
234 |
xpnt /= adj; |
235 |
ypnt /= adj; |
236 |
zpnt /= adj; |
237 |
xaxis( &xpnt, &ypnt, &zpnt ); |
238 |
if( xpnt > 0.F ) |
239 |
{ |
240 |
for( i3 = 1; i3 <= (natom + 4); i3++ ) |
241 |
atom[i3].x = -atom[i3].x; |
242 |
} |
243 |
/* **** after moving, check the # of h added to the carbon i */ |
244 |
if( nhadds == 1 ) |
245 |
{ |
246 |
/* **** there is only 1 h added to the carbon i */ |
247 |
for( i4 = 0; i4 < MAXIAT; i4++ ) |
248 |
{ |
249 |
if( atom[i].iat[i4] != 0 && atom[i].bo[i4] != 9 ) |
250 |
{ |
251 |
i5 = atom[i].iat[i4]; |
252 |
if( atom[i5].mmx_type == 5 || atom[i5].mmx_type == 20 ) |
253 |
{ |
254 |
/* i5 is the atom to attach */ |
255 |
if (adj < 2.0) // fix for OH in MMFF |
256 |
{ |
257 |
atom[i5].x = 0.70; |
258 |
atom[i5].y = 0.70; |
259 |
atom[i5].z = 0.00; |
260 |
} else |
261 |
{ |
262 |
atom[i5].x = 1.10F; |
263 |
if( atom[i5].mmx_type == 20 ) |
264 |
atom[i5].x = 0.6F; |
265 |
atom[i5].y = 0.00F; |
266 |
atom[i5].z = 0.00F; |
267 |
} |
268 |
if( atom[i].atomnum == 16 ) |
269 |
{ |
270 |
atom[i5].x = 0.F; |
271 |
atom[i5].y = 0.3796F; |
272 |
atom[i5].z = 1.043F; |
273 |
} |
274 |
} |
275 |
} |
276 |
} |
277 |
} |
278 |
/* **** there are 2 h added to the carbon i */ |
279 |
if( nhadds == 2 ) |
280 |
{ |
281 |
if( it == 2 || it == 3 || it == 29 || it == 30 || it == 9 || it == 41 || it == 46 |
282 |
|| it == 48 || it == 7 || it == 37 || it == 39 |
283 |
|| it == 42 || it == 38 || it == 40 || it == 44 || it == 66) |
284 |
{ |
285 |
for( i8 = 0; i8 < 10; i8++ ) |
286 |
{ |
287 |
if( iadj[i8] != 0 ) |
288 |
{ |
289 |
ii1 = iadj[i8]; |
290 |
for( ii2 = 0; ii2 < MAXIAT; ii2++ ) |
291 |
{ |
292 |
ii3 = 0; |
293 |
if( atom[ii1].iat[ii2] != 0 && atom[ii1].iat[ii2] != i && atom[ii1].bo[ii2] != 9 ) |
294 |
{ |
295 |
ii3 = atom[ii1].iat[ii2]; |
296 |
xpnt = atom[ii3].x; |
297 |
ypnt = atom[ii3].y; |
298 |
zpnt = atom[ii3].z; |
299 |
xyplan( &xpnt, &ypnt, &zpnt ); |
300 |
ij = 0; |
301 |
for( i9 = 0; i9 < MAXIAT; i9++ ) |
302 |
{ |
303 |
if( atom[i].iat[i9] != 0 && atom[i].bo[i9] != 9 ) |
304 |
{ |
305 |
iii = atom[i].iat[i9]; |
306 |
if( atom[iii].mmx_type == 5 || atom[iii].mmx_type == 20 ) |
307 |
{ |
308 |
ijk[ij] = iii; |
309 |
ij++; |
310 |
if( ij == 2 ) |
311 |
{ |
312 |
/* *** =ch2 *** */ |
313 |
atom[ijk[0]].x = 0.5582F; |
314 |
atom[ijk[1]].x = 0.5582F; |
315 |
atom[ijk[0]].y = 0.9478F; |
316 |
atom[ijk[1]].y = -0.9478F; |
317 |
atom[ijk[0]].z = 0.00F; |
318 |
atom[ijk[1]].z = 0.00F; |
319 |
if( it == 7 || it == 38 || it == 42 || it == 39 || it == 66) |
320 |
{ |
321 |
atom[ijk[0]].x = 0.30F; |
322 |
atom[ijk[1]].x = 0.30F; |
323 |
atom[ijk[0]].y = 0.52F; |
324 |
atom[ijk[1]].y = -0.52F; |
325 |
} |
326 |
if( atom[i].mmx_type == 37 ) |
327 |
{ |
328 |
atom[ijk[1]].x = 0.35F; |
329 |
atom[ijk[1]].y = -0.5F; |
330 |
} |
331 |
} |
332 |
} |
333 |
} |
334 |
} |
335 |
break; |
336 |
} |
337 |
if (ii3 == 0) // C=O with no lone pairs |
338 |
{ |
339 |
ij = 0; |
340 |
for( i9 = 0; i9 < MAXIAT; i9++ ) |
341 |
{ |
342 |
if( atom[i].iat[i9] != 0 && atom[i].bo[i9] != 9 ) |
343 |
{ |
344 |
iii = atom[i].iat[i9]; |
345 |
if( atom[iii].mmx_type == 5 || atom[iii].mmx_type == 20 ) |
346 |
{ |
347 |
ijk[ij] = iii; |
348 |
ij++; |
349 |
if( ij == 2 ) |
350 |
{ |
351 |
/* *** =ch2 *** */ |
352 |
atom[ijk[0]].x = 0.5582F; |
353 |
atom[ijk[1]].x = 0.5582F; |
354 |
atom[ijk[0]].y = 0.9478F; |
355 |
atom[ijk[1]].y = -0.9478F; |
356 |
atom[ijk[0]].z = 0.00F; |
357 |
atom[ijk[1]].z = 0.00F; |
358 |
if( it == 7 || it == 38 || it == 42 || it == 39 || it == 66) |
359 |
{ |
360 |
atom[ijk[0]].x = 0.30F; |
361 |
atom[ijk[1]].x = 0.30F; |
362 |
atom[ijk[0]].y = 0.52F; |
363 |
atom[ijk[1]].y = -0.52F; |
364 |
} |
365 |
if( atom[i].mmx_type == 37 ) |
366 |
{ |
367 |
atom[ijk[1]].x = 0.35F; |
368 |
atom[ijk[1]].y = -0.5F; |
369 |
} |
370 |
} |
371 |
} |
372 |
} |
373 |
} |
374 |
break; |
375 |
} |
376 |
} |
377 |
} |
378 |
} |
379 |
} else |
380 |
{ |
381 |
/* ** there are 2 atoms attached to carbon i, |
382 |
* one of which is "ni" and move these two atoms to xy-plane */ |
383 |
for( i6 = 0; i6 < 10; i6++ ) |
384 |
{ |
385 |
if( iadj[i6] != 0 ) |
386 |
{ |
387 |
ni = iadj[i6]; |
388 |
xpnt = atom[ni].x; |
389 |
ypnt = atom[ni].y; |
390 |
zpnt = atom[ni].z; |
391 |
xyplan( &xpnt, &ypnt, &zpnt ); |
392 |
} |
393 |
} |
394 |
ij = 0; |
395 |
for( i7 = 0; i7 < MAXIAT; i7++ ) |
396 |
{ |
397 |
if( atom[i].iat[i7] != 0 && atom[i].bo[i7] != 9 ) |
398 |
{ |
399 |
iii = atom[i].iat[i7]; |
400 |
if( atom[iii].mmx_type == 5 || atom[iii].mmx_type == 20 ) |
401 |
{ |
402 |
ijk[ij] = iii; |
403 |
ij++ ; |
404 |
if( ij == 2 ) |
405 |
{ |
406 |
/* *** -ch2- type 1 atom *** */ |
407 |
atom[ijk[0]].x = 0.6758F; |
408 |
atom[ijk[1]].x = 0.6758F; |
409 |
atom[ijk[0]].y = 0.00F; |
410 |
atom[ijk[1]].y = 0.00F; |
411 |
atom[ijk[0]].z = 0.8807F; |
412 |
atom[ijk[1]].z = -0.8807F; |
413 |
if( atom[i].mmx_type == 25 || atom[i].mmx_type == 46 || atom[i].mmx_type == 67) |
414 |
{ |
415 |
atom[ijk[0]].x = 0.3796F; |
416 |
atom[ijk[1]].x = 0.3796F; |
417 |
atom[ijk[0]].y = 1.043F; |
418 |
atom[ijk[1]].y = -0.5215F; |
419 |
atom[ijk[0]].z = 0.0000F; |
420 |
atom[ijk[1]].z = 0.9032F; |
421 |
} |
422 |
if( atom[i].mmx_type == 22 || atom[i].mmx_type == 30 ) |
423 |
{ |
424 |
/* *** -ch2- type 22 and 3o atoms *** */ |
425 |
atom[ijk[0]].x = 0.2F; |
426 |
atom[ijk[1]].x = 0.2F; |
427 |
atom[ijk[0]].y = 0.0F; |
428 |
atom[ijk[1]].y = 0.0F; |
429 |
atom[ijk[0]].z = 1.0F; |
430 |
atom[ijk[1]].z = -1.0F; |
431 |
} |
432 |
/* *** -o(lp)2- *** */ |
433 |
if (atom[i].mmx_type == 6 && (atom[iadj[0]].mmx_type == 3 || |
434 |
atom[iadj[0]].mmx_type == 2 || atom[iadj[0]].mmx_type == 30 || |
435 |
atom[iadj[0]].mmx_type == 37 || atom[iadj[0]].mmx_type == 40) ) |
436 |
{ |
437 |
atom[ijk[0]].x = 0.450F; |
438 |
atom[ijk[1]].x = 0.2488F; |
439 |
atom[ijk[0]].z = 0.0F; |
440 |
atom[ijk[1]].z = 0.0F; |
441 |
atom[ijk[0]].y = 0.9500F; |
442 |
atom[ijk[1]].y = -0.5460F; |
443 |
} |
444 |
else if( atom[i].mmx_type == 6 || atom[i].mmx_type == 15 || atom[i].mmx_type == |
445 |
42 || atom[i].mmx_type == 34 || atom[i].mmx_type == 35 || atom[i].mmx_type == 53 ) |
446 |
{ |
447 |
atom[ijk[0]].x = 0.2488F; |
448 |
atom[ijk[1]].x = 0.2488F; |
449 |
atom[ijk[0]].y = 0.0F; |
450 |
atom[ijk[1]].y = 0.0F; |
451 |
atom[ijk[0]].z = 0.5460F; |
452 |
atom[ijk[1]].z = -0.5460F; |
453 |
} |
454 |
} |
455 |
} |
456 |
} |
457 |
} |
458 |
} |
459 |
} |
460 |
/* end of nhadds = 2 |
461 |
* **** there are 3 h added to the carbon i *** */ |
462 |
if( nhadds == 3 ) |
463 |
{ |
464 |
for( j1 = 0; j1 < 10; j1++ ) |
465 |
{ |
466 |
if( iadj[j1] != 0 ) |
467 |
{ |
468 |
j2 = iadj[j1]; |
469 |
if( atom[j2].atomnum != 1 && atom[j2].mmx_type != 20 ) |
470 |
na = j2; |
471 |
if( atom[j2].mmx_type == 2 || atom[j2].mmx_type == 3 ) |
472 |
{ |
473 |
for( j3 = 0; j3 < MAXIAT; j3++ ) |
474 |
{ |
475 |
if( (atom[j2].iat[j3] != 0 && atom[j2].iat[j3] != i) && |
476 |
atom[j2].bo[j3] != 9 ) |
477 |
{ |
478 |
j4 = atom[j2].iat[j3]; |
479 |
if( atom[j4].mmx_type == 2 || (atom[j4].mmx_type == 7 && atom[j4].mmx_type == 38) ) |
480 |
{ |
481 |
xpnt = atom[j4].x; |
482 |
ypnt = atom[j4].y; |
483 |
zpnt = atom[j4].z; |
484 |
xyplan( &xpnt, &ypnt, &zpnt ); |
485 |
if( atom[j4].y < 0.F ) |
486 |
{ |
487 |
for( j5 = 1; j5 <= (natom + 4); j5++ ) |
488 |
atom[j5].y = -atom[j5].y; |
489 |
} |
490 |
} |
491 |
} |
492 |
} |
493 |
}else |
494 |
{ |
495 |
/* end of type 2 and type 3 adjacent* */ |
496 |
for( j8 = 0; j8 < MAXIAT; j8++ ) |
497 |
{ |
498 |
if( (atom[j2].iat[j8] != 0 && atom[j2].iat[j8] != i) |
499 |
&& atom[j2].bo[j8] != 9 ) |
500 |
{ |
501 |
j9 = atom[j2].iat[j8]; |
502 |
// if( atom[j9].mmx_type != 5 && atom[j9].mmx_type != 20 ) |
503 |
if(atom[j9].mmx_type != 20 ) |
504 |
{ |
505 |
xpnt = atom[j9].x; |
506 |
ypnt = atom[j9].y; |
507 |
zpnt = atom[j9].z; |
508 |
xyplan( &xpnt, &ypnt, &zpnt ); |
509 |
if( atom[j9].y > 0.F ) |
510 |
{ |
511 |
for( k1 = 1; k1 <= (natom + 4); k1++ ) |
512 |
atom[k1].y = -atom[k1].y; |
513 |
} |
514 |
} |
515 |
} |
516 |
} |
517 |
} |
518 |
/* do type 4 here */ |
519 |
ij = 0; |
520 |
for( k2 = 0; k2 < MAXIAT; k2++ ) |
521 |
{ |
522 |
if( !(atom[i].iat[k2] == 0 || atom[i].bo[k2] == 9) ) |
523 |
{ |
524 |
iii = atom[i].iat[k2]; |
525 |
if( !(atom[iii].mmx_type != 5 && atom[iii].mmx_type != 20) ) |
526 |
{ |
527 |
ijk[ij] = iii; |
528 |
ij++; |
529 |
if( ij == 3 ) |
530 |
{ |
531 |
dis = 0.F; |
532 |
/* *** -ch3 *** |
533 |
* *** add the hydrogens displaced along x enough |
534 |
* to give a normal length c-c bond *** */ |
535 |
atom[ijk[0]].x = 0.3796F + dis; |
536 |
atom[ijk[1]].x = 0.3796F + dis; |
537 |
atom[ijk[2]].x = 0.3796F + dis; |
538 |
atom[ijk[0]].y = 1.043F; |
539 |
atom[ijk[1]].y = -0.5215F; |
540 |
atom[ijk[2]].y = -0.5215F; |
541 |
atom[ijk[0]].z = 0.00F; |
542 |
atom[ijk[1]].z = 0.9032F; |
543 |
atom[ijk[2]].z = -0.9032F; |
544 |
} |
545 |
} |
546 |
} |
547 |
} |
548 |
} |
549 |
} |
550 |
} |
551 |
} |
552 |
} |
553 |
} |
554 |
} |
555 |
/* *** reorient the structure in starting orientation and |
556 |
* turn off atom markers *** */ |
557 |
revec(); |
558 |
return; |
559 |
} |
560 |
// ======================================================== |
561 |
void revec() |
562 |
{ |
563 |
int i; |
564 |
float xpnt, ypnt, zpnt; |
565 |
/* *** this routine used the unit vectors and origin saved |
566 |
* by 'subroutine getvec' to return the molecule to |
567 |
* its original orientation *** */ |
568 |
|
569 |
if( natom <= MAXATOM - 4 ) |
570 |
{ |
571 |
for( i = 1; i <= (natom + 4); i++ ) |
572 |
{ |
573 |
atom[i].x = atom[i].x - atom[natom + 4].x; |
574 |
atom[i].y = atom[i].y - atom[natom + 4].y; |
575 |
atom[i].z = atom[i].z - atom[natom + 4].z; |
576 |
} |
577 |
xpnt = atom[natom + 1].x; |
578 |
ypnt = atom[natom + 1].y; |
579 |
zpnt = atom[natom + 1].z; |
580 |
xaxis( &xpnt, &ypnt, &zpnt ); |
581 |
xpnt = atom[natom + 2].x; |
582 |
ypnt = atom[natom + 2].y; |
583 |
zpnt = atom[natom + 2].z; |
584 |
xyplan( &xpnt, &ypnt, &zpnt ); |
585 |
if( atom[natom + 1].x < 0. ) |
586 |
{ |
587 |
for( i = 1; i <= (natom + 4); i++ ) |
588 |
atom[i].x = -atom[i].x; |
589 |
} |
590 |
if( atom[natom + 2].y < 0. ) |
591 |
{ |
592 |
for( i = 1; i <= (natom + 4); i++ ) |
593 |
atom[i].y = -atom[i].y; |
594 |
} |
595 |
if( atom[natom + 3].z < 0. ) |
596 |
{ |
597 |
for( i = 1; i <= (natom + 4); i++ ) |
598 |
atom[i].z = -atom[i].z; |
599 |
} |
600 |
} |
601 |
return; |
602 |
} |
603 |
/* ------------------------ */ |
604 |
void xyplan(float *xr, float *yr, float *zr) |
605 |
{ |
606 |
int i; |
607 |
float cos3, denom, savex, sin3; |
608 |
|
609 |
/* this subroutine rotates the molecule to place atom (xr,yr,*zr) |
610 |
* in the x,y-plane. note0 y-coordinate will be +. */ |
611 |
|
612 |
denom = sqrt( *yr**yr + *zr**zr ); |
613 |
if( fabs( denom ) >= .00001F ) |
614 |
{ |
615 |
sin3 = *zr/denom; |
616 |
cos3 = *yr/denom; |
617 |
for( i = 1; i <= (natom + 4); i++ ) |
618 |
{ |
619 |
savex = atom[i].y; |
620 |
atom[i].y = atom[i].y*cos3 + atom[i].z*sin3; |
621 |
atom[i].z = atom[i].z*cos3 - savex*sin3; |
622 |
} |
623 |
savex = *yr; |
624 |
*yr = *yr*cos3 + *zr*sin3; |
625 |
*zr = *zr*cos3 - savex*sin3; |
626 |
} |
627 |
return; |
628 |
} |
629 |
/* ------------------------ */ |
630 |
void xaxis(float *xpnt,float *ypnt,float *zpnt) |
631 |
{ |
632 |
int i; |
633 |
float cos1, cos2, denom, save, sin1, sin2; |
634 |
|
635 |
/* this subroutine rotates the molecule about z-axis and then y-axis to |
636 |
* place point with coordinates (i*xpnt,iypnt,izpnt) along x-axis. */ |
637 |
|
638 |
denom = sqrt( *xpnt**xpnt + *ypnt**ypnt ); |
639 |
if( fabs( denom ) < .00001F ) |
640 |
{ |
641 |
sin1 = 0.F; |
642 |
cos1 = 1.F; |
643 |
}else |
644 |
{ |
645 |
sin1 = *ypnt/denom; |
646 |
cos1 = *xpnt/denom; |
647 |
for( i = 1; i <= (natom + 4); i++ ) |
648 |
{ |
649 |
save = atom[i].x; |
650 |
atom[i].x = atom[i].x*cos1 + atom[i].y*sin1; |
651 |
atom[i].y = atom[i].y*cos1 - save*sin1; |
652 |
} |
653 |
save = *xpnt; |
654 |
*xpnt = *xpnt*cos1 + *ypnt*sin1; |
655 |
*ypnt = *ypnt*cos1 - save*sin1; |
656 |
} |
657 |
denom = sqrt( *xpnt**xpnt + *zpnt**zpnt ); |
658 |
if( fabs( denom ) >= .00001F ) |
659 |
{ |
660 |
sin2 = *zpnt/denom; |
661 |
cos2 = *xpnt/denom; |
662 |
for( i = 1; i <= (natom + 4); i++ ) |
663 |
{ |
664 |
save = atom[i].x; |
665 |
atom[i].x = atom[i].x*cos2 + atom[i].z*sin2; |
666 |
atom[i].z = atom[i].z*cos2 - save*sin2; |
667 |
} |
668 |
save = *xpnt; |
669 |
*xpnt = *xpnt*cos2 + *zpnt*sin2; |
670 |
*zpnt = *zpnt*cos2 - save*sin2; |
671 |
} |
672 |
return; |
673 |
} |
674 |
/* ============================================== */ |
675 |
void getvec() |
676 |
{ |
677 |
|
678 |
/* *** this routine saves unit vectors and origin in the |
679 |
* n+1 - n+4 slots of x, y, and z *** */ |
680 |
|
681 |
if(natom > MAXATOM - 4 ) |
682 |
return; |
683 |
atom[natom + 1].x = 1.F; |
684 |
atom[natom + 1].y = 0.F; |
685 |
atom[natom + 1].z = 0.F; |
686 |
atom[natom + 2].x = 0.F; |
687 |
atom[natom + 2].y = 1.F; |
688 |
atom[natom + 2].z = 0.F; |
689 |
atom[natom + 3].x = 0.F; |
690 |
atom[natom + 3].y = 0.F; |
691 |
atom[natom + 3].z = 1.F; |
692 |
atom[natom + 4].x = 0.F; |
693 |
atom[natom + 4].y = 0.F; |
694 |
atom[natom + 4].z = 0.F; |
695 |
return; |
696 |
} |
697 |
/* ============================================== */ |
698 |
void hdel(int lptest) |
699 |
{ |
700 |
/* lptest = 0 to remove H and lone pair */ |
701 |
/* lptest = 1 to remove lone pairs only */ |
702 |
/* lptest = 2 to remove lp & reset various groups for amber */ |
703 |
int i, iatt, iatt1, iatyp, it, it1,ntemp; |
704 |
int j, jatom; |
705 |
/* *** delete hydrogens *** */ |
706 |
ntemp = natom; |
707 |
for( i = 1; i <= ntemp; i++ ) |
708 |
{ |
709 |
if( atom[i].mmx_type != 0 ) |
710 |
{ |
711 |
iatyp = atom[i].mmx_type; |
712 |
if( (lptest == 0 && (iatyp == 5 || iatyp == 20 )) || (lptest >= 1 && iatyp == 20) ) |
713 |
{ |
714 |
iatt = atom[i].iat[0]; |
715 |
iatt1 = atom[i].iat[1]; |
716 |
if( iatt != 0 || iatt1 != 0) |
717 |
{ |
718 |
it = atom[iatt].mmx_type; |
719 |
it1 = atom[iatt1].mmx_type; |
720 |
if (it >= 11 && it <= 14) // halogens |
721 |
goto L_10; |
722 |
if( (it < 300) && (it1 < 300)) // attached metals |
723 |
deleteatom(i); |
724 |
} |
725 |
} |
726 |
} |
727 |
L_10: |
728 |
continue; |
729 |
} |
730 |
reseq(); |
731 |
return; |
732 |
} |
733 |
/* -------------------------- */ |
734 |
void reseq() |
735 |
{ |
736 |
int i, ia, ib, iplus, ixx, j, jincr, k; |
737 |
|
738 |
/* ***this subroutine will repack the atom |
739 |
* ***connection table, making sure to fill in |
740 |
* **all "holes" caused by deleting atoms or bonds. |
741 |
* ***the necessary common blocks involve arrays |
742 |
* of the connection table (contb) */ |
743 |
|
744 |
for( i = 1; i < MAXATOM; i++ ) |
745 |
{ |
746 |
if( atom[i].type != 0 ) |
747 |
{ |
748 |
for( ia = 0; ia < (MAXIAT - 1); ia++ ) |
749 |
{ |
750 |
if( atom[i].iat[ia] == 0 ) |
751 |
{ |
752 |
/* *** for every zero entry in the array iat |
753 |
* for atom i, want to move down all remaining |
754 |
* entries of iat by one slot; after doing so |
755 |
* want to store that zero entry in the last |
756 |
* slot--iat *** */ |
757 |
for( ib = ia + 1; ib < MAXIAT; ib++ ) |
758 |
{ |
759 |
atom[i].iat[ib - 1] = atom[i].iat[ib]; |
760 |
atom[i].bo[ib - 1] = atom[i].bo[ib]; |
761 |
} |
762 |
atom[i].iat[MAXIAT - 1] = 0; |
763 |
} |
764 |
} |
765 |
} |
766 |
} |
767 |
|
768 |
/* *** now want to renumber atoms to leave |
769 |
* no "holes" caused by the deleted atoms *** */ |
770 |
jincr = 0; |
771 |
/* i points to first zero entry |
772 |
* kincr points to next real entry */ |
773 |
for( i = 1; i < (MAXATOM - 1); i++ ) |
774 |
{ |
775 |
iplus = 0; |
776 |
if( atom[i].type == 0 ) |
777 |
{ |
778 |
for( jincr = i; jincr < MAXATOM; jincr++ ) |
779 |
{ |
780 |
if( atom[jincr].type != 0 ) |
781 |
{ |
782 |
iplus = jincr; |
783 |
break; |
784 |
} |
785 |
} |
786 |
if (iplus == 0) |
787 |
break; |
788 |
for( j = 0; j < MAXIAT; j++ ) |
789 |
{ |
790 |
atom[i].iat[j] = atom[iplus].iat[j]; |
791 |
atom[i].bo[j] = atom[iplus].bo[j]; |
792 |
} |
793 |
|
794 |
atom[i].x = atom[iplus].x; |
795 |
atom[i].y = atom[iplus].y; |
796 |
atom[i].z = atom[iplus].z; |
797 |
atom[i].type = atom[iplus].type; |
798 |
atom[i].tclass = atom[iplus].tclass; |
799 |
atom[i].type = atom[iplus].type; |
800 |
atom[i].mmx_type = atom[iplus].mmx_type; |
801 |
atom[i].mm3_type = atom[iplus].mm3_type; |
802 |
atom[i].mmff_type = atom[iplus].mmff_type; |
803 |
atom[i].atomnum = atom[iplus].atomnum; |
804 |
atom[i].formal_charge = atom[iplus].formal_charge; |
805 |
atom[i].atomwt = atom[iplus].atomwt; |
806 |
atom[i].use = atom[iplus].use; |
807 |
strcpy(atom[i].name,atom[iplus].name); |
808 |
atom[i].charge = atom[iplus].charge; |
809 |
atom[i].flags = atom[iplus].flags; |
810 |
atom[i].radius = atom[iplus].radius; |
811 |
atom[i].energy = atom[iplus].energy; |
812 |
/* *** delete the atom whose statistics were just |
813 |
* stored in the i-th atom's connection table *** */ |
814 |
for( ixx = 0; ixx < MAXIAT; ixx++ ) |
815 |
{ |
816 |
atom[iplus].iat[ixx] = 0; |
817 |
atom[iplus].bo[ixx] = 0; |
818 |
} |
819 |
atom[iplus].x = 0.0F; |
820 |
atom[iplus].y = 0.0F; |
821 |
atom[iplus].z = 0.0F; |
822 |
atom[iplus].type = 0; |
823 |
atom[iplus].tclass = 0; |
824 |
atom[iplus].mmx_type = 0; |
825 |
atom[iplus].mm3_type = 0; |
826 |
atom[iplus].mmff_type = 0; |
827 |
atom[iplus].atomnum = 0; |
828 |
atom[iplus].formal_charge = 0; |
829 |
atom[iplus].atomwt = 0.0F; |
830 |
atom[iplus].use = 0; |
831 |
atom[iplus].charge = 0.0F; |
832 |
atom[iplus].flags = 0; |
833 |
atom[iplus].radius = 0.0F; |
834 |
atom[iplus].energy = 0.0F; |
835 |
strcpy(atom[iplus].name,""); |
836 |
|
837 |
for( j = 0; j < MAXIAT; j++ ) |
838 |
{ |
839 |
if( atom[i].iat[j] != 0 ) |
840 |
{ |
841 |
for( k = 0; k < MAXIAT; k++ ) |
842 |
{ |
843 |
if( atom[atom[i].iat[j]].iat[k] == iplus ) |
844 |
atom[atom[i].iat[j]].iat[k] = i; |
845 |
} |
846 |
} |
847 |
} |
848 |
} |
849 |
} |
850 |
/* ***now want to do the same resequencing for iat array */ |
851 |
|
852 |
/* recalc number of atoms */ |
853 |
natom = 0; |
854 |
for (i = 1; i < MAXATOM; i++) |
855 |
{ |
856 |
if (atom[i].type != 0) |
857 |
natom++; |
858 |
} |
859 |
return; |
860 |
} |