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 (13 years, 4 months ago) by wdelano
File size: 38228 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2     #include "pcwin.h"
3     #include "pcmod.h"
4     #include "field.h"
5     #include "atom_k.h"
6 gilbertke 63 #include "draw.h"
7 tjod 3
8     void hadd(void);
9     void hdel(int);
10     void type(void);
11     void nhadd(void);
12     void hcoord(void);
13 gilbertke 63 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 tjod 3
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 wdelano 99 // ========================
46 tjod 3 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 wdelano 99 set_atomdata(newatom,5,5,5);
99 tjod 3 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 wdelano 99 set_atomdata(newatom,20,0,0);
159 tjod 3 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 wdelano 58
837 tjod 3 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     }