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, 10 months ago) by wdelano
File size: 28010 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
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 double dihdrl(int i1,int i2,int i3,int i4);
8 void type_mmx(void);
9 void hdel(int);
10 void set_atomdata(int,int,int,int);
11
12 void set_atomtypes(int);
13 void type(void);
14 void quick_type(void);
15
16 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 }
52 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 } else if (jjbo < 4 && atom[i].mmx_type != 40 && atom[i].mmx_type != 4 && icount > 0)
88 {
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 set_atomdata(i,iatomt,0,mmfftype);
717 L_300:
718 continue;
719 }
720 }
721 return;
722 }
723