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