ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/pcmwin1.c
Revision: 96
Committed: Sat Jan 17 23:40:13 2009 UTC (10 years, 7 months ago) by wdelano
File size: 28985 byte(s)
Log Message:
using last pre-openmp source version
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,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 } 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 && atom[i].mmx_type != 4 && 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