ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/pcmwin2.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 6 months ago) by wdelano
File size: 38228 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
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 == 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 }