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

Line File contents
1 #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