ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/pcmwin1.c
Revision: 58
Committed: Mon Dec 1 06:44:59 2008 UTC (10 years, 10 months ago) by wdelano
Original Path: trunk/src/mengine/src/pcmwin1.c
File size: 28954 byte(s)
Log Message:
code merge 20081130
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     double dihdrl(int i1,int i2,int i3,int i4)
202     {
203     float a1, a2, ang, b1, b2, c1, c2, siign,
204     x1, x3, x4, y1, y3, y4, z1, z3, z4;
205     double dihdrl_v,r1,r2;
206    
207     /****the arguments of this function are the atom indices of atoms a-b-c-d
208     * which determine dihedral angle dihdrl.
209     * dihdrl returns the degree value of the dihedral('torsional')
210     * angle defined by atoms i1,i2,i3, &i4.
211     */
212    
213     dihdrl_v = 0.0;
214     if (i1 <= 0 || i2 <= 0 || i3 <= 0 || i4 <= 0)
215     return (dihdrl_v);
216     if (i1 > natom || i2 > natom || i3 > natom || i4 > natom)
217     return (dihdrl_v);
218    
219     a1 = atom[i2].x;
220     b1 = atom[i2].y;
221     c1 = atom[i2].z;
222     x1 = atom[i1].x - a1;
223     y1 = atom[i1].y - b1;
224     z1 = atom[i1].z - c1;
225     x3 = atom[i3].x - a1;
226     y3 = atom[i3].y - b1;
227     z3 = atom[i3].z - c1;
228     x4 = atom[i4].x - a1;
229     y4 = atom[i4].y - b1;
230     z4 = atom[i4].z - c1;
231     a1 = y1 * z3 - y3 * z1;
232     b1 = x3 * z1 - x1 * z3;
233     c1 = x1 * y3 - x3 * y1;
234     a2 = y4 * z3 - y3 * z4;
235     b2 = x3 * z4 - x4 * z3;
236     c2 = x4 * y3 - x3 * y4;
237     r1 = sqrt(a1 * a1 + b1 * b1 + c1 * c1);
238     r2 = sqrt(a2 * a2 + b2 * b2 + c2 * c2);
239     if (r1 != 0.0 && r2 != 0.0)
240     ang = (a1 * a2 + b1 * b2 + c1 * c2) /(r1*r2);
241     else
242     ang = 0.0;
243    
244     if (fabs(ang) > 1)
245     ang = fabs(ang) / ang;
246     ang = acos(ang);
247     dihdrl_v = 57.29578F * ang;
248     siign = x1 * a2 + y1 * b2 + z1 * c2;
249     if (siign < 0.0F)
250     dihdrl_v = -dihdrl_v;
251     return (dihdrl_v);
252     } /* end of function */
253     /* --------------------------------- */
254     /* ------------------------------------------------------------- */
255     void quick_type()
256     {
257     int i, iatik, iatkkk, iatomt, iatyp, ibt[MAXIAT], ibtk,
258     ibut, idc, itia, j, ja, jji, k, k1, k2, katm, ki,
259     kk, knt, l, lat, latm, lnt, k3, matm, jatm;
260     int mmfftype;
261    
262     /* *** loop over all atoms and determine type *** */
263     for (i = 1; i <= natom; i++)
264     {
265     if (atom[i].mmx_type != 0)
266     {
267     iatyp = atom[i].atomnum;
268     iatomt = atom[i].mmx_type;
269     /* *** load all attached bond types into ibt() *** except for
270     * coordinated bonds */
271     for (j = 0; j < MAXIAT; j++)
272     {
273     if (atom[i].bo[j] == 9 || atom[i].iat[j] == 0)
274     ibt[j] = 0;
275     else
276     ibt[j] = atom[i].bo[j];
277     }
278     /* *** now branch on the atom type *** */
279     if (iatyp == 6)
280     goto L_20;
281     if (iatyp == 1)
282     goto L_100;
283     if (iatyp == 8)
284     goto L_140;
285     if (iatyp == 7)
286     goto L_170;
287     if (iatyp == 16)
288     goto L_230;
289     if (iatyp == 34)
290     goto L_250;
291     if (iatyp == 15)
292     goto L_270;
293     if (iatyp == 13) // Aluminum
294     goto L_271;
295     if (iatyp == 17) // chlorine
296     goto L_272;
297     /* *** all other atoms have only one type *** */
298     goto L_300;
299    
300     /* *** carbon *** */
301     L_20:
302     iatomt = 1;
303     mmfftype = 1;
304     for (k = 0; k < MAXIAT; k++)
305     {
306     if (atom[i].iat[k] != 0 && atom[i].bo[k] >= 2 && atom[i].bo[k] != 9)
307     {
308     /* *** sp carbon *** */
309     if (atom[i].bo[k] == 3)
310     {
311     iatomt = 4;
312     mmfftype = 4;
313     goto L_290;
314     }
315     /* *** sp2 carbon *** */
316     ja = atom[i].iat[k];
317     iatomt = 2;
318     if (atom[ja].atomnum == 8 )
319     {
320     iatomt = 3;
321     mmfftype = 3;
322     /* check for ketenes and co2 */
323     kk = k + 1;
324     for (ki = kk; ki <= 10; ki++)
325     {
326     if (atom[i].bo[ki] == 2)
327     {
328     iatomt = 4;
329     mmfftype = 4;
330     goto L_290;
331     }
332     }
333     goto L_290;
334     }
335     /* check for allene */
336     kk = k + 1;
337     for (ki = kk; ki <= 8; ki++)
338     {
339     if (ibt[ki] == 2)
340     {
341     iatomt = 4;
342     mmfftype = 4;
343     goto L_290;
344     }
345     }
346     }
347     }
348     /* *** 3 and 4 -mem ring carbons *** */
349     if (iatomt != 1)
350     goto L_290;
351     for (k = 0; k < MAXIAT; k++)
352     {
353     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
354     {
355     katm = atom[i].iat[k];
356     for (k1 = 0; k1 < MAXIAT; k1++)
357     {
358     if (atom[katm].iat[k1] != 0 && atom[katm].iat[k1] != i && atom[katm].bo[k1] != 9)
359     {
360     latm = atom[katm].iat[k1];
361     for (k2 = 0; k2 < MAXIAT; k2++)
362     {
363     if (atom[latm].iat[k2] != 0 && atom[latm].iat[k2] != katm && atom[latm].bo[k2] != 9)
364     {
365     if (atom[latm].iat[k2] == i)
366     {
367     iatomt = 22;
368     mmfftype = 22;
369     ibut = 1;
370     goto L_290;
371     }
372     }
373     }
374     }
375     }
376     }
377     }
378     // now check for four membered rings
379     for (k = 0; k < MAXIAT; k++)
380     {
381     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
382     {
383     katm = atom[i].iat[k];
384     for (k1 = 0; k1 < MAXIAT; k1++)
385     {
386     if (atom[katm].iat[k1] != 0 && atom[katm].iat[k1] != i && atom[katm].bo[k1] != 9)
387     {
388     latm = atom[katm].iat[k1];
389     for (k2 = 0; k2 < MAXIAT; k2++)
390     {
391     if (atom[latm].iat[k2] != 0 && atom[latm].iat[k2] != katm && atom[latm].iat[k2] != i && atom[latm].bo[k2] != 9)
392     {
393     matm = atom[latm].iat[k2];
394     for (k3 = 0; k3 < MAXIAT; k3++)
395     {
396     if (atom[matm].iat[k3] != 0 && atom[matm].iat[k3] != latm && atom[matm].bo[k3] != 9)
397     {
398     if (atom[matm].iat[k3] == i)
399     {
400     iatomt = 56;
401     mmfftype = 20;
402     goto L_290;
403     }
404     }
405     }
406     }
407     }
408     }
409     }
410     }
411     }
412     goto L_290;
413    
414     /* *** hydrogen *** */
415     L_100:
416     iatomt = 5;
417     mmfftype = 5;
418     knt = 0;
419     for (k=0;k < MAXIAT; k++)
420     if (atom[i].iat[k] != 0)
421     knt++;
422     if (knt == 2)
423     {
424     iatomt = 70;
425     goto L_290;
426     }
427     for (k = 0; k < MAXIAT; k++)
428     {
429     if (atom[i].iat[k] != 0)
430     {
431     ja = atom[atom[i].iat[k]].mmx_type;
432     goto L_120;
433     }
434     }
435     goto L_290;
436     L_120:
437     if (((ja == 8 || ja == 9) || ja == 37) || ja == 55)
438     {
439     iatomt = 23;
440     mmfftype = 23;
441     } else if (ja == 41 || ja == 46)
442     {
443     iatomt = 24;
444     mmfftype = 36;
445     } else if (ja == 6 || ja == 53)
446     {
447     iatomt = 21;
448     mmfftype = 21;
449     for (j = 0; j < MAXIAT; j++) {
450     if (atom[atom[i].iat[k]].iat[j] == 0)
451     goto L_130;
452     if (atom[atom[i].iat[k]].iat[j] == i)
453     goto L_130;
454     if (atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 3)
455     {
456     iatomt = 24;
457     mmfftype = 24;
458     } else if (atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 2 ||
459     atom[atom[atom[i].iat[k]].iat[j]].mmx_type == 40)
460     {
461     iatomt = 28;
462     mmfftype = 28;
463     }
464     L_130:
465     ;
466     }
467     }
468     goto L_290;
469    
470     /* *** oxygen *** */
471     L_140:
472     if (iatomt == 42 || iatomt == 66)
473     {
474     mmfftype = 32;
475     goto L_290;
476     }
477     iatomt = 6;
478     mmfftype = 6;
479     for (k = 0; k < MAXIAT; k++)
480     {
481     if (atom[i].iat[k] != 0)
482     {
483     if (ibt[k] == 2)
484     {
485     iatomt = 7;
486     if (atom[atom[i].iat[0]].atomnum == 16)
487     mmfftype = 32;
488     goto L_290;
489     } else if (ibt[k] == 3)
490     {
491     for (lat = 0; lat < MAXIAT; lat++)
492     {
493     if (atom[i].iat[lat] != 0)
494     {
495     if (atom[atom[i].iat[lat]].mmx_type == 4)
496     {
497     iatomt = 46;
498     mmfftype = 49;
499     goto L_290;
500     }
501     }
502     }
503     }
504     }
505     }
506     goto L_290;
507    
508     /* *** nitrogen *** */
509     L_170:
510     iatomt = 8;
511     mmfftype = 8;
512     ibtk = 0;
513     /* determine maximum number of ligands attached to nitrogen */
514     jji = 0;
515     for (k = 0; k < MAXIAT; k++)
516     {
517     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
518     jji += 1;
519     }
520     for (k = 0; k < MAXIAT; k++)
521     {
522     if (ibt[k] != 9)
523     {
524     ibtk += ibt[k];
525     }
526     if (atom[i].iat[k] != 0)
527     {
528     if (atom[atom[i].iat[k]].mmx_type== 20 ||
529     (atom[atom[i].iat[k]].mmx_type >= 300) )
530     goto L_190;
531     }
532     }
533     /* ammonium nitrogen - no lone pairs */
534     if (ibtk == 4)
535     {
536     iatomt = 41;
537     mmfftype = 34;
538     goto L_290;
539     }
540     /* got a lone pair can not be ammonium */
541     L_190:
542     for (k = 0; k < MAXIAT; k++)
543     {
544     /* azo, imine */
545     if (ibt[k] == 2)
546     {
547     iatomt = 37;
548     mmfftype = 9;
549     goto L_290;
550     /* nitrile */
551     } else if (ibt[k] == 3)
552     {
553     iatomt = 10;
554     mmfftype = 42;
555     goto L_290;
556     }
557     }
558     /* count the number of double bonds to an adjacent Sulfur or
559     * Selenium */
560     idc = 0;
561     for (k = 0; k < MAXIAT; k++)
562     {
563     iatik = atom[i].iat[k];
564     if (iatik != 0 && atom[i].bo[k] != 9)
565     {
566     if (atom[iatik].atomnum == 16)
567     {
568     for (kk = 0; kk < MAXIAT; kk++)
569     {
570     iatkkk = atom[iatik].iat[kk];
571     if (iatkkk != 0 && atom[iatik].bo[kk] != 9)
572     {
573     if (atom[iatik].bo[kk] > 1)
574     idc += 1;
575     }
576     }
577     }
578     }
579     }
580     if (idc == 1)
581     {
582     iatomt = 9;
583     mmfftype = 9;
584     } else if (idc == 2)
585     {
586     iatomt = 8;
587     mmfftype = 9;
588     }
589     /* *** loop over all adjacent atoms to find enamines, amides */
590     for (k = 0; k < MAXIAT; k++)
591     {
592     if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9 )
593     {
594     itia = atom[atom[i].iat[k]].mmx_type;
595     if ( ((itia >= 300) || itia == 29 || itia == 30 || itia == 40 || itia == 26) && jji <= 3)
596     {
597     iatomt = 9;
598     mmfftype = 9;
599     goto L_290;
600     }
601     if (atom[atom[i].iat[k]].atomnum == 6) // adjacent carbon
602     {
603     jatm = atom[i].iat[k];
604     for (j=0; j < MAXIAT; j++)
605     {
606     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
607     {
608     iatomt = 9;
609     mmfftype = 10;
610     goto L_290;
611     }
612     }
613     }
614     if (atom[atom[i].iat[k]].atomnum == 16) // adjacent sulfur
615     {
616     jatm = atom[i].iat[k];
617     for (j=0; j < MAXIAT; j++)
618     {
619     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
620     {
621     iatomt = 9;
622     mmfftype = 9;
623     goto L_290;
624     }
625     }
626     }
627     if (atom[atom[i].iat[k]].atomnum == 7) // adjacent nitrogen
628     {
629     jatm = atom[i].iat[k];
630     for (j=0; j < MAXIAT; j++)
631     {
632     if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] >= 2)
633     {
634     iatomt = 9;
635     mmfftype = 9;
636     goto L_290;
637     } else if (atom[jatm].iat[j] != 0 && atom[jatm].iat[j] != i)
638     {
639     katm = atom[jatm].iat[j];
640     if (atom[katm].atomnum == 6)
641     {
642     for (l=0; l < MAXIAT; l++)
643     {
644     if (atom[katm].iat[l] != 0 && atom[katm].iat[l] != jatm && atom[katm].bo[l] == 2)
645     {
646     iatomt = 9;
647     mmfftype = 9;
648     goto L_290;
649     }
650     }
651     }
652     }
653     }
654     }
655     }
656     }
657     goto L_290;
658    
659     /* *** sulfur *** */
660     L_230:
661     if (iatomt == 16)
662     {
663     mmfftype = 16;
664     goto L_290;
665     }
666     iatomt = 15;
667     mmfftype = 15;
668     knt = 0;
669     lnt = 0;
670     for (k = 0; k < MAXIAT; k++) {
671     if (atom[i].iat[k] != 0) {
672     if (ibt[k] > 0 && atom[atom[i].iat[k]].mmx_type != 20)
673     {
674     lnt += 1;
675     }
676     if (ibt[k] > 1)
677     {
678     knt += 1;
679     if (atom[atom[i].iat[k]].mmx_type <= 4)
680     {
681     iatomt = 38;
682     mmfftype = 16;
683     goto L_290;
684     }
685     }
686     }
687     }
688     if (knt == 1)
689     {
690     iatomt = 17;
691     mmfftype = 17;
692     if (lnt == 1)
693     iatomt = 38;
694     }
695     if (knt >= 2)
696     {
697     iatomt = 18;
698     mmfftype = 18;
699     }
700     goto L_290;
701     /* ** selenium */
702     L_250:
703     iatomt = 34;
704     for (k = 0; k < MAXIAT; k++) {
705     if (atom[i].iat[k] != 0) {
706     if (ibt[k] == 2) {
707     iatomt = 39;
708     goto L_290;
709     }
710     }
711     }
712     goto L_290;
713     /* ** phosphorous */
714     L_270:
715     iatomt = 25;
716     mmfftype = 25;
717     knt = 0;
718     for (j = 0; j < MAXIAT; j++) {
719     if (atom[i].iat[j] != 0)
720     knt += 1;
721     }
722     if (knt >= 5)
723     iatomt = 47;
724     goto L_290;
725     /* aluminum */
726     L_271:
727     iatomt = 44;
728     knt = 0;
729     for (j = 0; j < MAXIAT; j++)
730     {
731     if (atom[i].iat[j] != 0)
732     knt += 1;
733     }
734     if (knt == 4)
735     iatomt = 58;
736     goto L_290;
737     /* Chlorine */
738     L_272:
739     iatomt = 12;
740     mmfftype = 12;
741     knt = 0;
742     for (j = 0; j < MAXIAT; j++)
743     {
744     if (atom[i].iat[j] != 0)
745     knt += 1;
746     }
747     if (knt == 2)
748     iatomt = 74;
749     goto L_290;
750     /* *** typing is now done, load value and return *** */
751     L_290:
752     // if (field.type == MMX && iatomt < 300 )
753     set_atomdata(i,iatomt,0,mmfftype,0,0);
754     L_300:
755     continue;
756     }
757     }
758     return;
759     }
760