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