ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/pcmwin2.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 7 months ago) by tjod
File size: 45824 byte(s)
Log Message:
test

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