ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/pcmwin2.c
Revision: 96
Committed: Sat Jan 17 23:40:13 2009 UTC (13 years, 4 months ago) by wdelano
File size: 38551 byte(s)
Log Message:
using last pre-openmp source version
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 == 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 wdelano 58
852 tjod 3 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     }