ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcmwin1.c
Revision: 63
Committed: Wed Dec 3 03:46:32 2008 UTC (10 years, 9 months ago) by gilbertke
File size: 28967 byte(s)
Log Message:
updated read_sdf type_mmx and first pass at using best practices
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 gilbertke 63 double dihdrl(int i1,int i2,int i3,int i4);
8 tjod 3 void type_mmx(void);
9     void hdel(int);
10     void set_atomdata(int,int,int,int,int,int);
11    
12 gilbertke 63 void set_atomtypes(int);
13     void type(void);
14     static void quick_type(void);
15    
16 tjod 3 void set_atomtypes(int newtype)
17     {
18     int i,badtype;
19     char hatext[256];
20    
21     if (newtype == MMX)
22     {
23     for (i=1; i <= natom; i++)
24     {
25     atom[i].type = atom[i].mmx_type;
26     if (atom[i].type < atom_k.natomtype)
27     atom[i].tclass = atom_k.tclass[atom[i].type];
28     else
29     atom[i].tclass = atom[i].type;
30     }
31     } else if (newtype == MM3)
32     {
33     for (i=1; i <= natom; i++)
34     {
35     atom[i].type = atom[i].mm3_type;
36     if (atom[i].type < atom_k.natomtype)
37     atom[i].tclass = atom_k.tclass[atom[i].type];
38     else
39     atom[i].tclass = atom[i].type;
40     }
41     } else if (newtype == MMFF94)
42     {
43     for (i=1; i <= natom; i++)
44     {
45     atom[i].type = atom[i].mmff_type;
46     if (atom[i].type < atom_k.natomtype)
47     atom[i].tclass = atom_k.tclass[atom[i].type];
48     else
49     atom[i].tclass = atom[i].type;
50     }
51     } else if (newtype == OPLSAA)
52     {
53     for (i=1; i <= natom; i++)
54     {
55     atom[i].type = atom[i].opls_type;
56     atom[i].tclass = atom_k.tclass[atom[i].type];
57     }
58     } else if (newtype == AMBER)
59     {
60     for (i=1; i <= natom; i++)
61     {
62     atom[i].type = atom[i].amber_type;
63     atom[i].tclass = atom_k.tclass[atom[i].type];
64     }
65     }
66     badtype = FALSE;
67    
68     for (i=1; i <= natom; i++)
69     if (atom[i].type == 0)
70     badtype = TRUE;
71     if (badtype == TRUE)
72     {
73     sprintf(hatext,"Error setting atom types in FField: %d \n",newtype);
74     message_alert(hatext,"ERROR");
75     }
76     }
77     void type()
78     {
79     int i, j, jjbo, icount,ncarbon, nhyd;
80    
81     icount = 0;
82     ncarbon = 0;
83     nhyd = 0;
84     for (i=1; i <= natom; i++)
85     {
86     if (atom[i].atomnum == 6)
87     {
88     jjbo = 0;
89     for (j=0; j < MAXIAT; j++)
90     {
91     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
92     jjbo += atom[i].bo[j];
93     if (atom[atom[i].iat[j]].atomnum == 1)
94     nhyd++;
95     }
96     if ( jjbo < 3 && atom[i].mmx_type == 40 && icount > 0 )
97     {
98     quick_type();
99     set_atomtypes(field.type);
100     return;
101     } else if (jjbo < 4 && atom[i].mmx_type != 40 && icount > 0)
102     {
103     quick_type();
104     set_atomtypes(field.type);
105     return;
106     }else
107     {
108     icount++;
109     ncarbon++;
110     }
111     if (icount > 15)
112     break;
113     } else if (atom[i].atomnum == 8)
114     {
115     jjbo = 0;
116     for (j=0; j < MAXIAT; j++)
117     {
118     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
119     jjbo += atom[i].bo[j];
120     }
121     if (jjbo < 1 && icount > 0)
122     {
123     quick_type();
124     set_atomtypes(field.type);
125     return;
126     } else if (jjbo == 1 && icount > 0 && atom[i].mmx_type == 6)
127     {
128     quick_type();
129     set_atomtypes(field.type);
130     return;
131     }else
132     icount++;
133     if (icount > 15)
134     break;
135     }
136    
137     }
138     if (ncarbon > 0 && nhyd == 0) // pathological cases where there are no carbons, or no carbons that normally have hydrogens
139     {
140     type_mmx();
141     set_atomtypes(field.type);
142     return;
143     }
144     // normal typing routines
145     if (field.type == MMX)
146     {
147     type_mmx();
148     set_atomtypes(MMX);
149     } else if (field.type == MM3)
150     {
151     for (i=1; i <= natom; i++)
152     {
153     if (atom[i].mmx_type == 20)
154     {
155     hdel(1);
156     break;
157     }
158     }
159     type_mmx();
160     set_atomtypes(MM3);
161     } else if (field.type == MMFF94)
162     {
163     for (i=1; i <= natom; i++)
164     {
165     if (atom[i].mmx_type == 20)
166     {
167     hdel(1);
168     break;
169     }
170     }
171     type_mmx();
172     set_atomtypes(MMFF94);
173     } else if (field.type == AMBER)
174     {
175     for (i=1; i <= natom; i++)
176     {
177     if (atom[i].mmx_type == 20)
178     {
179     hdel(1);
180     break;
181     }
182     }
183     type_mmx();
184     set_atomtypes(AMBER);
185     } else if (field.type == OPLSAA)
186     {
187     for (i=1; i <= natom; i++)
188     {
189     if (atom[i].mmx_type == 20)
190     {
191     hdel(1);
192     break;
193     }
194     }
195     type_mmx();
196     set_atomtypes(OPLSAA);
197     }
198     else
199     message_alert("No typing rules for this force field are implemented!","ERROR");
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     /* ------------------------------------------------------------- */
256     void quick_type()
257     {
258     int i, iatik, iatkkk, iatomt, iatyp, ibt[MAXIAT], ibtk,
259     ibut, idc, itia, j, ja, jji, k, k1, k2, katm, ki,
260     kk, knt, l, lat, latm, lnt, k3, matm, jatm;
261     int mmfftype;
262    
263     /* *** loop over all atoms and determine type *** */
264     for (i = 1; i <= natom; i++)
265     {
266     if (atom[i].mmx_type != 0)
267     {
268     iatyp = atom[i].atomnum;
269     iatomt = atom[i].mmx_type;
270     /* *** load all attached bond types into ibt() *** except for
271     * coordinated bonds */
272     for (j = 0; j < MAXIAT; j++)
273     {
274     if (atom[i].bo[j] == 9 || atom[i].iat[j] == 0)
275     ibt[j] = 0;
276     else
277     ibt[j] = atom[i].bo[j];
278     }
279     /* *** now branch on the atom type *** */
280     if (iatyp == 6)
281     goto L_20;
282     if (iatyp == 1)
283     goto L_100;
284     if (iatyp == 8)
285     goto L_140;
286     if (iatyp == 7)
287     goto L_170;
288     if (iatyp == 16)
289     goto L_230;
290     if (iatyp == 34)
291     goto L_250;
292     if (iatyp == 15)
293     goto L_270;
294     if (iatyp == 13) // Aluminum
295     goto L_271;
296     if (iatyp == 17) // chlorine
297     goto L_272;
298     /* *** all other atoms have only one type *** */
299     goto L_300;
300    
301     /* *** carbon *** */
302     L_20:
303     iatomt = 1;
304     mmfftype = 1;
305     for (k = 0; k < MAXIAT; k++)
306     {
307     if (atom[i].iat[k] != 0 && atom[i].bo[k] >= 2 && atom[i].bo[k] != 9)
308     {
309     /* *** sp carbon *** */
310     if (atom[i].bo[k] == 3)
311     {
312     iatomt = 4;
313     mmfftype = 4;
314     goto L_290;
315     }
316     /* *** sp2 carbon *** */
317     ja = atom[i].iat[k];
318     iatomt = 2;
319     if (atom[ja].atomnum == 8 )
320     {
321     iatomt = 3;
322     mmfftype = 3;
323     /* check for ketenes and co2 */
324     kk = k + 1;
325     for (ki = kk; ki <= 10; ki++)
326     {
327     if (atom[i].bo[ki] == 2)
328     {
329     iatomt = 4;
330     mmfftype = 4;
331     goto L_290;
332     }
333     }
334     goto L_290;
335     }
336     /* check for allene */
337     kk = k + 1;
338     for (ki = kk; ki <= 8; ki++)
339     {
340     if (ibt[ki] == 2)
341     {
342     iatomt = 4;
343     mmfftype = 4;
344     goto L_290;
345     }
346     }
347     }
348     }
349     /* *** 3 and 4 -mem ring carbons *** */
350     if (iatomt != 1)
351     goto L_290;
352     for (k = 0; k < MAXIAT; k++)
353     {
354     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
355     {
356     katm = atom[i].iat[k];
357     for (k1 = 0; k1 < MAXIAT; k1++)
358     {
359     if (atom[katm].iat[k1] != 0 && atom[katm].iat[k1] != i && atom[katm].bo[k1] != 9)
360     {
361     latm = atom[katm].iat[k1];
362     for (k2 = 0; k2 < MAXIAT; k2++)
363     {
364     if (atom[latm].iat[k2] != 0 && atom[latm].iat[k2] != katm && atom[latm].bo[k2] != 9)
365     {
366     if (atom[latm].iat[k2] == i)
367     {
368     iatomt = 22;
369     mmfftype = 22;
370     ibut = 1;
371     goto L_290;
372     }
373     }
374     }
375     }
376     }
377     }
378     }
379     // now check for four membered rings
380     for (k = 0; k < MAXIAT; k++)
381     {
382     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
383     {
384     katm = atom[i].iat[k];
385     for (k1 = 0; k1 < MAXIAT; k1++)
386     {
387     if (atom[katm].iat[k1] != 0 && atom[katm].iat[k1] != i && atom[katm].bo[k1] != 9)
388     {
389     latm = atom[katm].iat[k1];
390     for (k2 = 0; k2 < MAXIAT; k2++)
391     {
392     if (atom[latm].iat[k2] != 0 && atom[latm].iat[k2] != katm && atom[latm].iat[k2] != i && atom[latm].bo[k2] != 9)
393     {
394     matm = atom[latm].iat[k2];
395     for (k3 = 0; k3 < MAXIAT; k3++)
396     {
397     if (atom[matm].iat[k3] != 0 && atom[matm].iat[k3] != latm && atom[matm].bo[k3] != 9)
398     {
399     if (atom[matm].iat[k3] == i)
400     {
401     iatomt = 56;
402     mmfftype = 20;
403     goto L_290;
404     }
405     }
406     }
407     }
408     }
409     }
410     }
411     }
412     }
413     goto L_290;
414    
415     /* *** hydrogen *** */
416     L_100:
417     iatomt = 5;
418     mmfftype = 5;
419     knt = 0;
420     for (k=0;k < MAXIAT; k++)
421     if (atom[i].iat[k] != 0)
422     knt++;
423     if (knt == 2)
424     {
425     iatomt = 70;
426     goto L_290;
427     }
428     for (k = 0; k < MAXIAT; k++)
429     {
430     if (atom[i].iat[k] != 0)
431     {
432     ja = atom[atom[i].iat[k]].mmx_type;
433     goto L_120;
434     }
435     }
436     goto L_290;
437     L_120:
438     if (((ja == 8 || ja == 9) || ja == 37) || ja == 55)
439     {
440     iatomt = 23;
441     mmfftype = 23;
442     } else if (ja == 41 || ja == 46)
443     {
444     iatomt = 24;
445     mmfftype = 36;
446     } else if (ja == 6 || ja == 53)
447     {
448     iatomt = 21;
449     mmfftype = 21;
450     for (j = 0; j < MAXIAT; j++) {
451     if (atom[atom[i].iat[k]].iat[j] == 0)
452     goto L_130;
453     if (atom[atom[i].iat[k]].iat[j] == i)
454     goto L_130;
455     if (atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 3)
456     {
457     iatomt = 24;
458     mmfftype = 24;
459     } else if (atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 2 ||
460     atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 40)
461     {
462     iatomt = 28;
463     mmfftype = 28;
464     }
465     L_130:
466     ;
467     }
468     }
469     goto L_290;
470    
471     /* *** oxygen *** */
472     L_140:
473     if (iatomt == 42 || iatomt == 66)
474     {
475     mmfftype = 32;
476     goto L_290;
477     }
478     iatomt = 6;
479     mmfftype = 6;
480     for (k = 0; k < MAXIAT; k++)
481     {
482     if (atom[i].iat[k] != 0)
483     {
484     if (ibt[k] == 2)
485     {
486     iatomt = 7;
487     if (atom[atom[i].iat[0]].atomnum == 16)
488     mmfftype = 32;
489     goto L_290;
490     } else if (ibt[k] == 3)
491     {
492     for (lat = 0; lat < MAXIAT; lat++)
493     {
494     if (atom[i].iat[lat] != 0)
495     {
496     if (atom[atom[i].iat[lat]].mmx_type == 4)
497     {
498     iatomt = 46;
499     mmfftype = 49;
500     goto L_290;
501     }
502     }
503     }
504     }
505     }
506     }
507     goto L_290;
508    
509     /* *** nitrogen *** */
510     L_170:
511     iatomt = 8;
512     mmfftype = 8;
513     ibtk = 0;
514     /* determine maximum number of ligands attached to nitrogen */
515     jji = 0;
516     for (k = 0; k < MAXIAT; k++)
517     {
518     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
519     jji += 1;
520     }
521     for (k = 0; k < MAXIAT; k++)
522     {
523     if (ibt[k] != 9)
524     {
525     ibtk += ibt[k];
526     }
527     if (atom[i].iat[k] != 0)
528     {
529     if (atom[atom[i].iat[k]].mmx_type== 20 ||
530     (atom[atom[i].iat[k]].mmx_type >= 300) )
531     goto L_190;
532     }
533     }
534     /* ammonium nitrogen - no lone pairs */
535     if (ibtk == 4)
536     {
537     iatomt = 41;
538     mmfftype = 34;
539     goto L_290;
540     }
541     /* got a lone pair can not be ammonium */
542     L_190:
543     for (k = 0; k < MAXIAT; k++)
544     {
545     /* azo, imine */
546     if (ibt[k] == 2)
547     {
548     iatomt = 37;
549     mmfftype = 9;
550     goto L_290;
551     /* nitrile */
552     } else if (ibt[k] == 3)
553     {
554     iatomt = 10;
555     mmfftype = 42;
556     goto L_290;
557     }
558     }
559     /* count the number of double bonds to an adjacent Sulfur or
560     * Selenium */
561     idc = 0;
562     for (k = 0; k < MAXIAT; k++)
563     {
564     iatik = atom[i].iat[k];
565     if (iatik != 0 && atom[i].bo[k] != 9)
566     {
567     if (atom[iatik].atomnum == 16)
568     {
569     for (kk = 0; kk < MAXIAT; kk++)
570     {
571     iatkkk = atom[iatik].iat[kk];
572     if (iatkkk != 0 && atom[iatik].bo[kk] != 9)
573     {
574     if (atom[iatik].bo[kk] > 1)
575     idc += 1;
576     }
577     }
578     }
579     }
580     }
581     if (idc == 1)
582     {
583     iatomt = 9;
584     mmfftype = 9;
585     } else if (idc == 2)
586     {
587     iatomt = 8;
588     mmfftype = 9;
589     }
590     /* *** loop over all adjacent atoms to find enamines, amides */
591     for (k = 0; k < MAXIAT; k++)
592     {
593     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9 )
594     {
595     itia = atom[atom[i].iat[k]].mmx_type;
596     if ( ((itia >= 300) || itia == 29 || itia == 30 || itia == 40 || itia == 26) && jji <= 3)
597     {
598     iatomt = 9;
599     mmfftype = 9;
600     goto L_290;
601     }
602     if (atom[atom[i].iat[k]].atomnum == 6) // adjacent carbon
603     {
604     jatm = atom[i].iat[k];
605     for (j=0; j < MAXIAT; j++)
606     {
607     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
608     {
609     iatomt = 9;
610     mmfftype = 10;
611     goto L_290;
612     }
613     }
614     }
615     if (atom[atom[i].iat[k]].atomnum == 16) // adjacent sulfur
616     {
617     jatm = atom[i].iat[k];
618     for (j=0; j < MAXIAT; j++)
619     {
620     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
621     {
622     iatomt = 9;
623     mmfftype = 9;
624     goto L_290;
625     }
626     }
627     }
628     if (atom[atom[i].iat[k]].atomnum == 7) // adjacent nitrogen
629     {
630     jatm = atom[i].iat[k];
631     for (j=0; j < MAXIAT; j++)
632     {
633     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
634     {
635     iatomt = 9;
636     mmfftype = 9;
637     goto L_290;
638     } else if (atom[jatm].iat[j] != 0 && atom[jatm].iat[j] != i)
639     {
640     katm = atom[jatm].iat[j];
641     if (atom[katm].atomnum == 6)
642     {
643     for (l=0; l < MAXIAT; l++)
644     {
645     if (atom[katm].iat[l] != 0 && atom[katm].iat[l] != jatm && atom[katm].bo[l] == 2)
646     {
647     iatomt = 9;
648     mmfftype = 9;
649     goto L_290;
650     }
651     }
652     }
653     }
654     }
655     }
656     }
657     }
658     goto L_290;
659    
660     /* *** sulfur *** */
661     L_230:
662     if (iatomt == 16)
663     {
664     mmfftype = 16;
665     goto L_290;
666     }
667     iatomt = 15;
668     mmfftype = 15;
669     knt = 0;
670     lnt = 0;
671     for (k = 0; k < MAXIAT; k++) {
672     if (atom[i].iat[k] != 0) {
673     if (ibt[k] > 0 && atom[atom[i].iat[k]].mmx_type != 20)
674     {
675     lnt += 1;
676     }
677     if (ibt[k] > 1)
678     {
679     knt += 1;
680     if (atom[atom[i].iat[k]].mmx_type <= 4)
681     {
682     iatomt = 38;
683     mmfftype = 16;
684     goto L_290;
685     }
686     }
687     }
688     }
689     if (knt == 1)
690     {
691     iatomt = 17;
692     mmfftype = 17;
693     if (lnt == 1)
694     iatomt = 38;
695     }
696     if (knt >= 2)
697     {
698     iatomt = 18;
699     mmfftype = 18;
700     }
701     goto L_290;
702     /* ** selenium */
703     L_250:
704     iatomt = 34;
705     for (k = 0; k < MAXIAT; k++) {
706     if (atom[i].iat[k] != 0) {
707     if (ibt[k] == 2) {
708     iatomt = 39;
709     goto L_290;
710     }
711     }
712     }
713     goto L_290;
714     /* ** phosphorous */
715     L_270:
716     iatomt = 25;
717     mmfftype = 25;
718     knt = 0;
719     for (j = 0; j < MAXIAT; j++) {
720     if (atom[i].iat[j] != 0)
721     knt += 1;
722     }
723     if (knt >= 5)
724     iatomt = 47;
725     goto L_290;
726     /* aluminum */
727     L_271:
728     iatomt = 44;
729     knt = 0;
730     for (j = 0; j < MAXIAT; j++)
731     {
732     if (atom[i].iat[j] != 0)
733     knt += 1;
734     }
735     if (knt == 4)
736     iatomt = 58;
737     goto L_290;
738     /* Chlorine */
739     L_272:
740     iatomt = 12;
741     mmfftype = 12;
742     knt = 0;
743     for (j = 0; j < MAXIAT; j++)
744     {
745     if (atom[i].iat[j] != 0)
746     knt += 1;
747     }
748     if (knt == 2)
749     iatomt = 74;
750     goto L_290;
751     /* *** typing is now done, load value and return *** */
752     L_290:
753     // if (field.type == MMX && iatomt < 300 )
754     set_atomdata(i,iatomt,0,mmfftype,0,0);
755     L_300:
756     continue;
757     }
758     }
759     return;
760     }
761