ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/pcmwin1.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (11 years, 4 months ago) by tjod
Original Path: trunk/src/mengine/src/pcmwin1.c
File size: 37627 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

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 type_mmx(void);
8     void type_mm3(void);
9     void hdel(int);
10     void quick_type(void);
11     void set_atomdata(int,int,int,int,int,int);
12     void message_alert(char *, char *);
13     void set_atomtypes(int);
14    
15     void set_atomtypes(int newtype)
16     {
17     int i,badtype;
18     char hatext[256];
19    
20     if (newtype == MMX)
21     {
22     for (i=1; i <= natom; i++)
23     {
24     atom[i].type = atom[i].mmx_type;
25     if (atom[i].type < atom_k.natomtype)
26     atom[i].tclass = atom_k.tclass[atom[i].type];
27     else
28     atom[i].tclass = atom[i].type;
29     }
30     } else if (newtype == MM3)
31     {
32     for (i=1; i <= natom; i++)
33     {
34     atom[i].type = atom[i].mm3_type;
35     if (atom[i].type < atom_k.natomtype)
36     atom[i].tclass = atom_k.tclass[atom[i].type];
37     else
38     atom[i].tclass = atom[i].type;
39     }
40     } else if (newtype == MMFF94)
41     {
42     for (i=1; i <= natom; i++)
43     {
44     atom[i].type = atom[i].mmff_type;
45     if (atom[i].type < atom_k.natomtype)
46     atom[i].tclass = atom_k.tclass[atom[i].type];
47     else
48     atom[i].tclass = atom[i].type;
49     }
50     } else if (newtype == OPLSAA)
51     {
52     for (i=1; i <= natom; i++)
53     {
54     atom[i].type = atom[i].opls_type;
55     atom[i].tclass = atom_k.tclass[atom[i].type];
56     }
57     } else if (newtype == AMBER)
58     {
59     for (i=1; i <= natom; i++)
60     {
61     atom[i].type = atom[i].amber_type;
62     atom[i].tclass = atom_k.tclass[atom[i].type];
63     }
64     }
65     badtype = FALSE;
66    
67     for (i=1; i <= natom; i++)
68     if (atom[i].type == 0)
69     badtype = TRUE;
70     if (badtype == TRUE)
71     {
72     sprintf(hatext,"Error setting atom types in FField: %d \n",newtype);
73     message_alert(hatext,"ERROR");
74     }
75     }
76     void type()
77     {
78     int i, j, jjbo, icount,ncarbon, nhyd;
79    
80     icount = 0;
81     ncarbon = 0;
82     nhyd = 0;
83     for (i=1; i <= natom; i++)
84     {
85     if (atom[i].atomnum == 6)
86     {
87     jjbo = 0;
88     for (j=0; j < MAXIAT; j++)
89     {
90     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
91     jjbo += atom[i].bo[j];
92     if (atom[atom[i].iat[j]].atomnum == 1)
93     nhyd++;
94     }
95     if ( jjbo < 3 && atom[i].mmx_type == 40 && icount > 0 )
96     {
97     quick_type();
98     set_atomtypes(field.type);
99     return;
100     } else if (jjbo < 4 && atom[i].mmx_type != 40 && icount > 0)
101     {
102     quick_type();
103     set_atomtypes(field.type);
104     return;
105     }else
106     {
107     icount++;
108     ncarbon++;
109     }
110     if (icount > 15)
111     break;
112     } else if (atom[i].atomnum == 8)
113     {
114     jjbo = 0;
115     for (j=0; j < MAXIAT; j++)
116     {
117     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
118     jjbo += atom[i].bo[j];
119     }
120     if (jjbo < 1 && icount > 0)
121     {
122     quick_type();
123     set_atomtypes(field.type);
124     return;
125     } else if (jjbo == 1 && icount > 0 && atom[i].mmx_type == 6)
126     {
127     quick_type();
128     set_atomtypes(field.type);
129     return;
130     }else
131     icount++;
132     if (icount > 15)
133     break;
134     }
135    
136     }
137     if (ncarbon > 0 && nhyd == 0) // pathological cases where there are no carbons, or no carbons that normally have hydrogens
138     {
139     type_mmx();
140     set_atomtypes(field.type);
141     return;
142     }
143     // normal typing routines
144     if (field.type == MMX)
145     {
146     type_mmx();
147     set_atomtypes(MMX);
148     } else if (field.type == MM3)
149     {
150     for (i=1; i <= natom; i++)
151     {
152     if (atom[i].mmx_type == 20)
153     {
154     hdel(1);
155     break;
156     }
157     }
158     type_mmx();
159     set_atomtypes(MM3);
160     } else if (field.type == MMFF94)
161     {
162     for (i=1; i <= natom; i++)
163     {
164     if (atom[i].mmx_type == 20)
165     {
166     hdel(1);
167     break;
168     }
169     }
170     type_mmx();
171     set_atomtypes(MMFF94);
172     } else if (field.type == AMBER)
173     {
174     for (i=1; i <= natom; i++)
175     {
176     if (atom[i].mmx_type == 20)
177     {
178     hdel(1);
179     break;
180     }
181     }
182     type_mmx();
183     set_atomtypes(AMBER);
184     } else if (field.type == OPLSAA)
185     {
186     for (i=1; i <= natom; i++)
187     {
188     if (atom[i].mmx_type == 20)
189     {
190     hdel(1);
191     break;
192     }
193     }
194     type_mmx();
195     set_atomtypes(OPLSAA);
196     }
197     else
198     message_alert("No typing rules for this force field are implemented!","ERROR");
199     }
200     /* ------------------------ */
201    
202     double dihdrl(int i1,int i2,int i3,int i4)
203     {
204     float a1, a2, ang, b1, b2, c1, c2, siign,
205     x1, x3, x4, y1, y3, y4, z1, z3, z4;
206     double dihdrl_v,r1,r2;
207    
208     /****the arguments of this function are the atom indices of atoms a-b-c-d
209     * which determine dihedral angle dihdrl.
210     * dihdrl returns the degree value of the dihedral('torsional')
211     * angle defined by atoms i1,i2,i3, &i4.
212     */
213    
214     dihdrl_v = 0.0;
215     if (i1 <= 0 || i2 <= 0 || i3 <= 0 || i4 <= 0)
216     return (dihdrl_v);
217     if (i1 > natom || i2 > natom || i3 > natom || i4 > natom)
218     return (dihdrl_v);
219    
220     a1 = atom[i2].x;
221     b1 = atom[i2].y;
222     c1 = atom[i2].z;
223     x1 = atom[i1].x - a1;
224     y1 = atom[i1].y - b1;
225     z1 = atom[i1].z - c1;
226     x3 = atom[i3].x - a1;
227     y3 = atom[i3].y - b1;
228     z3 = atom[i3].z - c1;
229     x4 = atom[i4].x - a1;
230     y4 = atom[i4].y - b1;
231     z4 = atom[i4].z - c1;
232     a1 = y1 * z3 - y3 * z1;
233     b1 = x3 * z1 - x1 * z3;
234     c1 = x1 * y3 - x3 * y1;
235     a2 = y4 * z3 - y3 * z4;
236     b2 = x3 * z4 - x4 * z3;
237     c2 = x4 * y3 - x3 * y4;
238     r1 = sqrt(a1 * a1 + b1 * b1 + c1 * c1);
239     r2 = sqrt(a2 * a2 + b2 * b2 + c2 * c2);
240     if (r1 != 0.0 && r2 != 0.0)
241     ang = (a1 * a2 + b1 * b2 + c1 * c2) /(r1*r2);
242     else
243     ang = 0.0;
244    
245     if (fabs(ang) > 1)
246     ang = fabs(ang) / ang;
247     ang = acos(ang);
248     dihdrl_v = 57.29578F * ang;
249     siign = x1 * a2 + y1 * b2 + z1 * c2;
250     if (siign < 0.0F)
251     dihdrl_v = -dihdrl_v;
252     return (dihdrl_v);
253     } /* end of function */
254     /* --------------------------------- */
255     void avgleg()
256     {
257     int i, iatyp, ib1, ib2, j;
258     int hasSulf;
259     float an, xx, xyzdis, yy, zz;
260    
261     /**** if there is no bond, average=1.0 *** */
262    
263     an = 0.0F;
264     xyzdis = 0.0F;
265     hasSulf = FALSE;
266     for( i = 1; i <= natom; i++ )
267     {
268     if( atom[i].type && atom[i].type != 20 )
269     {
270     for( j = 0; j < MAXIAT; j++ )
271     {
272     if( atom[i].iat[j] && atom[i].bo[j] != 9 )
273     {
274     /* count bonds only once - don't count bonds to hydrogen */
275     if( atom[i].iat[j] >= i )
276     {
277     ib1 = i;
278     iatyp = atom[i].type;
279     if (atom[i].atomnum == 16) hasSulf = TRUE;
280     if( iatyp != 5 && iatyp != 21 && iatyp != 23 && iatyp != 24 && iatyp != 28 && iatyp != 36 && iatyp != 15 && iatyp != 38 )
281     {
282     ib2 = atom[i].iat[j];
283     iatyp = atom[atom[i].iat[j]].type;
284     if( iatyp != 20 )
285     {
286     if( iatyp != 5 && iatyp != 21 && iatyp != 23 && iatyp != 24 && iatyp != 28 && iatyp != 36 && iatyp != 15 && iatyp != 38)
287     {
288     an += 1.0F;
289     xx = atom[ib1].x - atom[ib2].x;
290     yy = atom[ib1].y - atom[ib2].y;
291     zz = atom[ib1].z - atom[ib2].z;
292     xyzdis += sqrt( xx*xx + yy*yy + zz*zz );
293     }
294     }
295     }
296     }
297     }
298     }
299     }
300     }
301     /* if an = 0 then there are no heavy atom bonds
302     * need to redo calculations based on c-h bonds */
303     if( an != 0.0F )
304     {
305     if (hasSulf)
306     flags.avleg = 1.73F/(xyzdis/an); /* assume normal c-c distance */
307     else
308     flags.avleg = 1.53F/(xyzdis/an); /* assume normal c-c distance */
309     }
310     else if( natom > 1 )
311     {
312     xyzdis = 0.0F;
313     an = 0.0F;
314     for( i = 1; i <= natom; i++ )
315     {
316     ib1 = i;
317     for( j = 0; j < MAXIAT; j++ )
318     {
319     if( atom[i].iat[j] && atom[i].bo[j] != 9 )
320     {
321     ib2 = atom[i].iat[j];
322     if( ib2 > ib1 )
323     {
324     an += 1.0F;
325     xx = atom[ib1].x - atom[ib2].x;
326     yy = atom[ib1].y - atom[ib2].y;
327     zz = atom[ib1].z - atom[ib2].z;
328     xyzdis += sqrt( xx*xx + yy*yy + zz*zz );
329     }
330     }
331     }
332     }
333     if (an > 0.0)
334     flags.avleg = 1.10F/(xyzdis/an); /* no c-c bonds assume c-h bond dist */
335     else
336     flags.avleg = 1.20F; /* default when all else goes wrong */
337     }
338     else
339     {
340     flags.avleg = 1.20F; /* default when all else goes wrong */
341     }
342     return;
343     }
344     /* -------------------------- */
345     void center(int istruct)
346     {
347     // istruct = -2 do all
348     // istruct = 1 do main structrure
349     // istruct = number do that substructure
350    
351     float centerx, centery, centerz;
352     int i;
353    
354     centerx = 0.0F;
355     centery = 0.0F;
356     centerz = 0.0F;
357    
358     if (istruct == -2)
359     {
360     for (i=1; i <= natom; i++)
361     {
362     centerx += atom[i].x/natom;
363     centery += atom[i].y/natom;
364     centerz += atom[i].z/natom;
365     }
366     for (i=1; i <= natom; i++)
367     {
368     atom[i].x -= centerx;
369     atom[i].y -= centery;
370     atom[i].z -= centerz;
371     }
372     } else if (istruct == -1) // main structure only
373     {
374     for (i=1; i <= natom; i++)
375     {
376     if(atom[i].substr[0] == 0)
377     {
378     centerx += atom[i].x/natom;
379     centery += atom[i].y/natom;
380     centerz += atom[i].z/natom;
381     }
382     }
383     for (i=1; i <= natom; i++)
384     {
385     if(atom[i].substr[0] == 0)
386     {
387     atom[i].x -= centerx;
388     atom[i].y -= centery;
389     atom[i].z -= centerz;
390     }
391     }
392     } else // substructure
393     {
394     for (i=1; i <= natom; i++)
395     {
396     if (atom[i].substr[0] & (1L << istruct) )
397     {
398     centerx += atom[i].x/natom;
399     centery += atom[i].y/natom;
400     centerz += atom[i].z/natom;
401     }
402     }
403     for (i=1; i <= natom; i++)
404     {
405     if (atom[i].substr[0] & (1L << istruct) )
406     {
407     atom[i].x -= centerx;
408     atom[i].y -= centery;
409     atom[i].z -= centerz;
410     }
411     }
412     }
413     }
414     /* ------------------------ */
415    
416     double ooplane(int ia,int ib,int ic,int id)
417     {
418     double angle,dot,cosine;
419     double xia,yia,zia;
420     double xib,yib,zib;
421     double xic,yic,zic;
422     double xid,yid,zid;
423     double xab,yab,zab,rab2,rab;
424     double xbc,ybc,zbc,rbc2,rbc;
425     double xbd,ybd,zbd,rbd2,rbd;
426     double sine;
427    
428     /****the arguments of this function are the atom indices of atoms a-b-c-d
429     * which determine dihedral angle dihdrl.
430     * dihdrl returns the degree value of the dihedral('torsional')
431     * angle defined by atoms i1,i2,i3, &i4.
432     */
433    
434     xia = atom[ia].x;
435     yia = atom[ia].y;
436     zia = atom[ia].z;
437     xib = atom[ib].x;
438     yib = atom[ib].y;
439     zib = atom[ib].z;
440     xic = atom[ic].x;
441     yic = atom[ic].y;
442     zic = atom[ic].z;
443     xid = atom[id].x;
444     yid = atom[id].y;
445     zid = atom[id].z;
446    
447     xab = xia - xib;
448     yab = yia - yib;
449     zab = zia - zib;
450     rab2 = (xab*xab+yab*yab+zab*zab);
451     rab = sqrt(rab2);
452    
453     xbc = xic - xib;
454     ybc = yic - yib;
455     zbc = zic - zib;
456     rbc2 = xbc*xbc+ybc*ybc+zbc*zbc;
457     rbc = sqrt(rbc2);
458    
459     if (rab2 != 0.0 && rbc2 != 0.0)
460     {
461     dot = xab*xbc + yab*ybc + zab*zbc;
462     cosine = dot /(rab*rbc);
463     if (cosine < -1.0)
464     cosine = -1.0;
465     if (cosine > 1.0)
466     cosine = 1.0;
467     sine = sqrt(1.00-cosine*cosine);
468    
469     xbd = atom[id].x - atom[ib].x;
470     ybd = atom[id].y - atom[ib].y;
471     zbd = atom[id].z - atom[ib].z;
472     rbd2 = xbd*xbd+ybd*ybd+zbd*zbd;
473     rbd = sqrt(rbd2);
474    
475     cosine = ( xbd*(yab*zbc-zab*ybc) + ybd*(zab*xbc-xab*zbc) + zbd*(xab*ybc-yab*xbc))/((rab*rbc*rbd)*sine);
476    
477     angle = radian * asin(cosine);
478     }
479     return (angle);
480     }
481     /* ------------------------------------------------------------- */
482     void quick_type()
483     {
484     int i, iatik, iatkkk, iatomt, iatyp, ibt[MAXIAT], ibtk,
485     ibut, idc, itia, j, ja, jji, k, k1, k2, katm, ki,
486     kk, knt, l, lat, latm, lnt, k3, matm, jatm;
487     int mmfftype;
488    
489     /* *** loop over all atoms and determine type *** */
490     for (i = 1; i <= natom; i++)
491     {
492     if (atom[i].mmx_type != 0)
493     {
494     iatyp = atom[i].atomnum;
495     iatomt = atom[i].mmx_type;
496     /* *** load all attached bond types into ibt() *** except for
497     * coordinated bonds */
498     for (j = 0; j < MAXIAT; j++)
499     {
500     if (atom[i].bo[j] == 9 || atom[i].iat[j] == 0)
501     ibt[j] = 0;
502     else
503     ibt[j] = atom[i].bo[j];
504     }
505     /* *** now branch on the atom type *** */
506     if (iatyp == 6)
507     goto L_20;
508     if (iatyp == 1)
509     goto L_100;
510     if (iatyp == 8)
511     goto L_140;
512     if (iatyp == 7)
513     goto L_170;
514     if (iatyp == 16)
515     goto L_230;
516     if (iatyp == 34)
517     goto L_250;
518     if (iatyp == 15)
519     goto L_270;
520     if (iatyp == 13) // Aluminum
521     goto L_271;
522     if (iatyp == 17) // chlorine
523     goto L_272;
524     /* *** all other atoms have only one type *** */
525     goto L_300;
526    
527     /* *** carbon *** */
528     L_20:
529     iatomt = 1;
530     mmfftype = 1;
531     for (k = 0; k < MAXIAT; k++)
532     {
533     if (atom[i].iat[k] != 0 && atom[i].bo[k] >= 2 && atom[i].bo[k] != 9)
534     {
535     /* *** sp carbon *** */
536     if (atom[i].bo[k] == 3)
537     {
538     iatomt = 4;
539     mmfftype = 4;
540     goto L_290;
541     }
542     /* *** sp2 carbon *** */
543     ja = atom[i].iat[k];
544     iatomt = 2;
545     if (atom[ja].atomnum == 8 )
546     {
547     iatomt = 3;
548     mmfftype = 3;
549     /* check for ketenes and co2 */
550     kk = k + 1;
551     for (ki = kk; ki <= 10; ki++)
552     {
553     if (atom[i].bo[ki] == 2)
554     {
555     iatomt = 4;
556     mmfftype = 4;
557     goto L_290;
558     }
559     }
560     goto L_290;
561     }
562     /* check for allene */
563     kk = k + 1;
564     for (ki = kk; ki <= 8; ki++)
565     {
566     if (ibt[ki] == 2)
567     {
568     iatomt = 4;
569     mmfftype = 4;
570     goto L_290;
571     }
572     }
573     }
574     }
575     /* *** 3 and 4 -mem ring carbons *** */
576     if (iatomt != 1)
577     goto L_290;
578     for (k = 0; k < MAXIAT; k++)
579     {
580     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
581     {
582     katm = atom[i].iat[k];
583     for (k1 = 0; k1 < MAXIAT; k1++)
584     {
585     if (atom[katm].iat[k1] != 0 && atom[katm].iat[k1] != i && atom[katm].bo[k1] != 9)
586     {
587     latm = atom[katm].iat[k1];
588     for (k2 = 0; k2 < MAXIAT; k2++)
589     {
590     if (atom[latm].iat[k2] != 0 && atom[latm].iat[k2] != katm && atom[latm].bo[k2] != 9)
591     {
592     if (atom[latm].iat[k2] == i)
593     {
594     iatomt = 22;
595     mmfftype = 22;
596     ibut = 1;
597     goto L_290;
598     }
599     }
600     }
601     }
602     }
603     }
604     }
605     // now check for four membered rings
606     for (k = 0; k < MAXIAT; k++)
607     {
608     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
609     {
610     katm = atom[i].iat[k];
611     for (k1 = 0; k1 < MAXIAT; k1++)
612     {
613     if (atom[katm].iat[k1] != 0 && atom[katm].iat[k1] != i && atom[katm].bo[k1] != 9)
614     {
615     latm = atom[katm].iat[k1];
616     for (k2 = 0; k2 < MAXIAT; k2++)
617     {
618     if (atom[latm].iat[k2] != 0 && atom[latm].iat[k2] != katm && atom[latm].iat[k2] != i && atom[latm].bo[k2] != 9)
619     {
620     matm = atom[latm].iat[k2];
621     for (k3 = 0; k3 < MAXIAT; k3++)
622     {
623     if (atom[matm].iat[k3] != 0 && atom[matm].iat[k3] != latm && atom[matm].bo[k3] != 9)
624     {
625     if (atom[matm].iat[k3] == i)
626     {
627     iatomt = 56;
628     mmfftype = 20;
629     goto L_290;
630     }
631     }
632     }
633     }
634     }
635     }
636     }
637     }
638     }
639     goto L_290;
640    
641     /* *** hydrogen *** */
642     L_100:
643     iatomt = 5;
644     mmfftype = 5;
645     knt = 0;
646     for (k=0;k < MAXIAT; k++)
647     if (atom[i].iat[k] != 0)
648     knt++;
649     if (knt == 2)
650     {
651     iatomt = 70;
652     goto L_290;
653     }
654     for (k = 0; k < MAXIAT; k++)
655     {
656     if (atom[i].iat[k] != 0)
657     {
658     ja = atom[atom[i].iat[k]].mmx_type;
659     goto L_120;
660     }
661     }
662     goto L_290;
663     L_120:
664     if (((ja == 8 || ja == 9) || ja == 37) || ja == 55)
665     {
666     iatomt = 23;
667     mmfftype = 23;
668     } else if (ja == 41 || ja == 46)
669     {
670     iatomt = 24;
671     mmfftype = 36;
672     } else if (ja == 6 || ja == 53)
673     {
674     iatomt = 21;
675     mmfftype = 21;
676     for (j = 0; j < MAXIAT; j++) {
677     if (atom[atom[i].iat[k]].iat[j] == 0)
678     goto L_130;
679     if (atom[atom[i].iat[k]].iat[j] == i)
680     goto L_130;
681     if (atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 3)
682     {
683     iatomt = 24;
684     mmfftype = 24;
685     } else if (atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 2 ||
686     atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 40)
687     {
688     iatomt = 28;
689     mmfftype = 28;
690     }
691     L_130:
692     ;
693     }
694     }
695     goto L_290;
696    
697     /* *** oxygen *** */
698     L_140:
699     if (iatomt == 42 || iatomt == 66)
700     {
701     mmfftype = 32;
702     goto L_290;
703     }
704     iatomt = 6;
705     mmfftype = 6;
706     for (k = 0; k < MAXIAT; k++)
707     {
708     if (atom[i].iat[k] != 0)
709     {
710     if (ibt[k] == 2)
711     {
712     iatomt = 7;
713     if (atom[atom[i].iat[0]].atomnum == 16)
714     mmfftype = 32;
715     goto L_290;
716     } else if (ibt[k] == 3)
717     {
718     for (lat = 0; lat < MAXIAT; lat++)
719     {
720     if (atom[i].iat[lat] != 0)
721     {
722     if (atom[atom[i].iat[lat]].mmx_type == 4)
723     {
724     iatomt = 46;
725     mmfftype = 49;
726     goto L_290;
727     }
728     }
729     }
730     }
731     }
732     }
733     goto L_290;
734    
735     /* *** nitrogen *** */
736     L_170:
737     iatomt = 8;
738     mmfftype = 8;
739     ibtk = 0;
740     /* determine maximum number of ligands attached to nitrogen */
741     jji = 0;
742     for (k = 0; k < MAXIAT; k++)
743     {
744     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
745     jji += 1;
746     }
747     for (k = 0; k < MAXIAT; k++)
748     {
749     if (ibt[k] != 9)
750     {
751     ibtk += ibt[k];
752     }
753     if (atom[i].iat[k] != 0)
754     {
755     if (atom[atom[i].iat[k]].mmx_type== 20 ||
756     (atom[atom[i].iat[k]].mmx_type >= 300) )
757     goto L_190;
758     }
759     }
760     /* ammonium nitrogen - no lone pairs */
761     if (ibtk == 4)
762     {
763     iatomt = 41;
764     mmfftype = 34;
765     goto L_290;
766     }
767     /* got a lone pair can not be ammonium */
768     L_190:
769     for (k = 0; k < MAXIAT; k++)
770     {
771     /* azo, imine */
772     if (ibt[k] == 2)
773     {
774     iatomt = 37;
775     mmfftype = 9;
776     goto L_290;
777     /* nitrile */
778     } else if (ibt[k] == 3)
779     {
780     iatomt = 10;
781     mmfftype = 42;
782     goto L_290;
783     }
784     }
785     /* count the number of double bonds to an adjacent Sulfur or
786     * Selenium */
787     idc = 0;
788     for (k = 0; k < MAXIAT; k++)
789     {
790     iatik = atom[i].iat[k];
791     if (iatik != 0 && atom[i].bo[k] != 9)
792     {
793     if (atom[iatik].atomnum == 16)
794     {
795     for (kk = 0; kk < MAXIAT; kk++)
796     {
797     iatkkk = atom[iatik].iat[kk];
798     if (iatkkk != 0 && atom[iatik].bo[kk] != 9)
799     {
800     if (atom[iatik].bo[kk] > 1)
801     idc += 1;
802     }
803     }
804     }
805     }
806     }
807     if (idc == 1)
808     {
809     iatomt = 9;
810     mmfftype = 9;
811     } else if (idc == 2)
812     {
813     iatomt = 8;
814     mmfftype = 9;
815     }
816     /* *** loop over all adjacent atoms to find enamines, amides */
817     for (k = 0; k < MAXIAT; k++)
818     {
819     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9 )
820     {
821     itia = atom[atom[i].iat[k]].mmx_type;
822     if ( ((itia >= 300) || itia == 29 || itia == 30 || itia == 40 || itia == 26) && jji <= 3)
823     {
824     iatomt = 9;
825     mmfftype = 9;
826     goto L_290;
827     }
828     if (atom[atom[i].iat[k]].atomnum == 6) // adjacent carbon
829     {
830     jatm = atom[i].iat[k];
831     for (j=0; j < MAXIAT; j++)
832     {
833     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
834     {
835     iatomt = 9;
836     mmfftype = 10;
837     goto L_290;
838     }
839     }
840     }
841     if (atom[atom[i].iat[k]].atomnum == 16) // adjacent sulfur
842     {
843     jatm = atom[i].iat[k];
844     for (j=0; j < MAXIAT; j++)
845     {
846     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
847     {
848     iatomt = 9;
849     mmfftype = 9;
850     goto L_290;
851     }
852     }
853     }
854     if (atom[atom[i].iat[k]].atomnum == 7) // adjacent nitrogen
855     {
856     jatm = atom[i].iat[k];
857     for (j=0; j < MAXIAT; j++)
858     {
859     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
860     {
861     iatomt = 9;
862     mmfftype = 9;
863     goto L_290;
864     } else if (atom[jatm].iat[j] != 0 && atom[jatm].iat[j] != i)
865     {
866     katm = atom[jatm].iat[j];
867     if (atom[katm].atomnum == 6)
868     {
869     for (l=0; l < MAXIAT; l++)
870     {
871     if (atom[katm].iat[l] != 0 && atom[katm].iat[l] != jatm && atom[katm].bo[l] == 2)
872     {
873     iatomt = 9;
874     mmfftype = 9;
875     goto L_290;
876     }
877     }
878     }
879     }
880     }
881     }
882     }
883     }
884     goto L_290;
885    
886     /* *** sulfur *** */
887     L_230:
888     if (iatomt == 16)
889     {
890     mmfftype = 16;
891     goto L_290;
892     }
893     iatomt = 15;
894     mmfftype = 15;
895     knt = 0;
896     lnt = 0;
897     for (k = 0; k < MAXIAT; k++) {
898     if (atom[i].iat[k] != 0) {
899     if (ibt[k] > 0 && atom[atom[i].iat[k]].mmx_type != 20)
900     {
901     lnt += 1;
902     }
903     if (ibt[k] > 1)
904     {
905     knt += 1;
906     if (atom[atom[i].iat[k]].mmx_type <= 4)
907     {
908     iatomt = 38;
909     mmfftype = 16;
910     goto L_290;
911     }
912     }
913     }
914     }
915     if (knt == 1)
916     {
917     iatomt = 17;
918     mmfftype = 17;
919     if (lnt == 1)
920     iatomt = 38;
921     }
922     if (knt >= 2)
923     {
924     iatomt = 18;
925     mmfftype = 18;
926     }
927     goto L_290;
928     /* ** selenium */
929     L_250:
930     iatomt = 34;
931     for (k = 0; k < MAXIAT; k++) {
932     if (atom[i].iat[k] != 0) {
933     if (ibt[k] == 2) {
934     iatomt = 39;
935     goto L_290;
936     }
937     }
938     }
939     goto L_290;
940     /* ** phosphorous */
941     L_270:
942     iatomt = 25;
943     mmfftype = 25;
944     knt = 0;
945     for (j = 0; j < MAXIAT; j++) {
946     if (atom[i].iat[j] != 0)
947     knt += 1;
948     }
949     if (knt >= 5)
950     iatomt = 47;
951     goto L_290;
952     /* aluminum */
953     L_271:
954     iatomt = 44;
955     knt = 0;
956     for (j = 0; j < MAXIAT; j++)
957     {
958     if (atom[i].iat[j] != 0)
959     knt += 1;
960     }
961     if (knt == 4)
962     iatomt = 58;
963     goto L_290;
964     /* Chlorine */
965     L_272:
966     iatomt = 12;
967     mmfftype = 12;
968     knt = 0;
969     for (j = 0; j < MAXIAT; j++)
970     {
971     if (atom[i].iat[j] != 0)
972     knt += 1;
973     }
974     if (knt == 2)
975     iatomt = 74;
976     goto L_290;
977     /* *** typing is now done, load value and return *** */
978     L_290:
979     // if (field.type == MMX && iatomt < 300 )
980     set_atomdata(i,iatomt,0,mmfftype,0,0);
981     L_300:
982     continue;
983     }
984     }
985     return;
986     }
987