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