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, 8 months ago) by wdelano
Original Path: trunk/src/mengine/src/pcmwin1.c
File size: 28954 byte(s)
Log Message:
code merge 20081130
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 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