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