ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/kangle.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 11 months ago) by tjod
File size: 24249 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2     #include "pcwin.h"
3     #include "pcmod.h"
4    
5     #include "angles.h"
6     #include "rings.h"
7     #include "field.h"
8     #include "bonds_ff.h"
9     #include "atom_k.h"
10    
11     int is_ring43(int, int, int);
12     int is_ring53(int, int, int);
13     int is_ring63(int, int, int);
14     int is_delocalbond(int, int);
15     int is_TS(int, int, int);
16     int isbond(int,int);
17     int find_bond(int,int);
18     void transtheta(int *, float *, float *, char *, char *, char *, char *);
19     void ts_four(int *, float *, float *,char *, char *,char *,char *);
20     long ipack(int ,int ,int ,int);
21     long kang(int,int,int);
22     void numeral(int,char *,int);
23     float get_general_angle(int);
24    
25    
26     EXTERN struct t_ts_bondorder {
27     float fbnd[15];
28     } ts_bondorder;
29    
30     EXTERN struct t_angk1 {
31     int use_ang3, use_ang4, use_ang5;
32     int nang, nang3, nang4, nang5;
33     int ndel, ndel3, ndel4;
34     char ktype[MAXANGCONST][10],ktype3[MAXANG3CONST][10],ktype4[MAXANG4CONST][10],ktype5[MAXANG5CONST][10];
35     char kdel[MAXANGDEL][10],kdel3[MAXANG3DEL][10],kdel4[MAXANG4DEL][10];
36     float con[MAXANGCONST], ang[MAXANGCONST][3];
37     float con3[MAXANG3CONST], ang3[MAXANG3CONST][3];
38     float con4[MAXANG4CONST], ang4[MAXANG4CONST][3];
39     float con5[MAXANG5CONST], ang5[MAXANG5CONST][3];
40     float condel[MAXANGDEL], angdel[MAXANGDEL][3];
41     float condel3[MAXANG3DEL], angdel3[MAXANG3DEL][3];
42     float condel4[MAXANG4DEL], angdel4[MAXANG4DEL][3];
43     } angk1;
44    
45     EXTERN struct t_angkp1 {
46     int npiang;
47     char kpa[MAXPIANGCONST][10];
48     float pacon[MAXPIANGCONST], panat[MAXPIANGCONST][3], piobp[MAXPIANGCONST];
49     } angkp1;
50    
51     EXTERN struct t_angf {
52     int use_angf, nfang;
53     char kftype[MAXANGCONST][10];
54     float fcon[MAXANGCONST],fc0[MAXANGCONST],fc1[MAXANGCONST],fc2[MAXANGCONST];
55     float fc3[MAXANGCONST],fc4[MAXANGCONST],fc5[MAXANGCONST],fc6[MAXANGCONST];
56     } angf;
57    
58     EXTERN int Missing_constants;
59    
60     void kangle()
61     {
62     long int pi_mask, aromatic_mask;
63     int i, j, k, izero, ierr, nh, nRc, nligand, jji;
64     int ia, ib, ic, ie;
65     int iaa, ibb, icc;
66     int cl_a, cl_b, cl_c;
67     int nbo1, nb1, nb2;
68     char pa[4],pb[4],pc[4],pt[10], pt1[10], pt2[10],pt3[10], pt4[10], pt5[10], pt6[10], pt7[10];
69     char zero[4];
70     float bkan, at;
71    
72     pi_mask = (1L << PI_MASK);
73     aromatic_mask = (1 << AROMATIC_MASK);
74    
75     if( angles.nang != 0 )
76     {
77     for( i = 0; i < angles.nang; i++ )
78     {
79     ia = angles.i13[i][0];
80     ib = angles.i13[i][1];
81     ic = angles.i13[i][2];
82     iaa = atom[ia].type;
83     ibb = atom[ib].type;
84     icc = atom[ic].type;
85     cl_a = atom[ia].tclass;
86     cl_b = atom[ib].tclass;
87     cl_c = atom[ic].tclass;
88    
89     /* get bond index each bond for MMFF94 */
90     nb1 = find_bond(ia,ib);
91     nb2 = find_bond(ib,ic);
92    
93     /* special metal cases */
94     if( iaa >= 300)
95     cl_a = 300;
96    
97     if( ibb >= 300)
98     cl_b = 300;
99    
100     if( icc >= 300)
101     cl_c = 300;
102    
103     /* special MMX cases */
104     if ( field.type == MMX)
105     {
106     if( iaa == 40 )
107     iaa = 2;
108    
109     if( ibb == 40 )
110     ibb = 2;
111    
112     if( icc == 40 )
113     icc = 2;
114    
115     }
116     numeral(iaa,pa,3);
117     numeral(ibb,pb,3);
118     numeral(icc,pc,3);
119     if( iaa <= icc )
120     {
121     strcpy(pt,pa);
122     strcat(pt,pb);
123     strcat(pt,pc);
124     }
125     if( iaa > icc )
126     {
127     strcpy(pt,pc);
128     strcat(pt,pb);
129     strcat(pt,pa);
130     }
131     izero = 0;
132     strcpy(zero," 0");
133     strcpy(pt1,pa); strcat(pt1,pb); strcat(pt1,zero);
134     strcpy(pt2,pc); strcat(pt2,pb); strcat(pt2,zero);
135     strcpy(pt3,zero); strcat(pt3,pb); strcat(pt3,pc);
136     strcpy(pt4,zero); strcat(pt4,pb); strcat(pt4,pa);
137     strcpy(pt5,zero); strcat(pt5,pb); strcat(pt5,zero);
138    
139     numeral(cl_a,pa,3);
140     numeral(cl_b,pb,3);
141     numeral(cl_c,pc,3);
142     if (cl_a < cl_c)
143     {
144     strcpy(pt6,pa); strcat(pt6,pb); strcat(pt6,pc);
145     }else
146     {
147     strcpy(pt6,pc); strcat(pt6,pb); strcat(pt6,pa);
148     }
149    
150     if (field.type == MMFF94)
151     {
152     cl_a = atom_k.tclass1[atom[ia].type];
153     cl_b = atom_k.tclass1[atom[ib].type];
154     cl_c = atom_k.tclass1[atom[ic].type];
155     numeral(cl_a,pa,3);
156     numeral(cl_b,pb,3);
157     numeral(cl_c,pc,3);
158     if (cl_a < cl_c)
159     {
160     strcpy(pt7,pa); strcat(pt7,pb); strcat(pt7,pc);
161     }else
162     {
163     strcpy(pt7,pc); strcat(pt7,pb); strcat(pt7,pa);
164     }
165     }
166     angles.acon[i] = 0.0;
167     angles.anat[i] = 0.0;
168     angles.angtype[i] = HARMONIC;
169    
170     ierr = FALSE;
171     /* check TS constants */
172     /* check delocalized angles in MMFF94 cyclopropanes */
173     if (field.type == MMFF94 && angk1.ndel3 > 0)
174     {
175     if (isbond(ia,ic) )
176     {
177     nbo1 = bonds_ff.index[nb1] + bonds_ff.index[nb2];
178     if (nbo1 == 1 || nbo1 == 2)
179     {
180     for(j=0; j < angk1.ndel3; j++)
181     {
182     if (strcmp(angk1.kdel3[j],pt) == 0)
183     {
184     angles.acon[i] = angk1.condel3[j];
185     angles.anat[i] = angk1.angdel3[j][0];
186     if (nbo1 == 1)
187     angles.index[i] = 5;
188     else if (nbo1 == 2)
189     angles.index[i] = 6;
190     ierr = TRUE;
191     break;
192     }
193     }
194     }
195    
196     if (ierr == TRUE)
197     goto L_10;
198     }
199     }
200     /* check three membered rings */
201     if ( (rings.nring3 > 0) && (angk1.nang3 > 0) )
202     {
203     if (isbond(ia,ic))
204     {
205     /* search three membered ring parameters */
206     for (j = 0; j < angk1.nang3; j++)
207     {
208     if ( (strcmp(angk1.ktype3[j],pt) == 0) || (strcmp(angk1.ktype3[j],pt1) == 0) ||
209     (strcmp(angk1.ktype3[j],pt2) == 0) || (strcmp(angk1.ktype3[j],pt3) == 0) ||
210     (strcmp(angk1.ktype3[j],pt4) == 0) || (strcmp(angk1.ktype3[j],pt5) == 0))
211     {
212     if ( angk1.ang3[j][1] == 0.0 && angk1.ang3[j][2] == 0)
213     {
214     angles.acon[i] = angk1.con3[j];
215     angles.anat[i] = angk1.ang3[j][0];
216     angles.index[i] = 3;
217     ierr = TRUE;
218     break;
219     } else
220     {
221     nh = 0;
222     for (k=0; k < MAXIAT; k++)
223     {
224     ie = atom[ib].iat[k];
225     if ( (atom[ie].type == 5 || atom[ie].type == 36) && (ie != ia) && (ie != ic))
226     nh++;
227     }
228     angles.acon[i] = angk1.con3[j];
229     angles.anat[i] = angk1.ang3[j][nh];
230     angles.index[i] = 3;
231     ierr = TRUE;
232     break;
233     }
234     }
235     }
236     if (ierr == TRUE)
237     goto L_10;
238     else
239     {
240     Missing_constants = TRUE;
241     }
242     }
243     }
244     /* check delocalized angles in MMFF94 cyclobutanes */
245     if (field.type == MMFF94 && angk1.ndel4 > 0)
246     {
247     nRc = is_ring43(ia, ib, ic);
248     if (nRc == TRUE)
249     {
250     nbo1 = bonds_ff.index[nb1] + bonds_ff.index[nb2];
251     if (nbo1 == 1 || nbo1 == 2)
252     {
253     for(j=0; j < angk1.ndel4; j++)
254     {
255     if (strcmp(angk1.kdel4[j],pt) == 0)
256     {
257     angles.acon[i] = angk1.condel4[j];
258     angles.anat[i] = angk1.angdel4[j][0];
259     if (nbo1 == 1)
260     angles.index[i] = 7;
261     else if (nbo1 == 2)
262     angles.index[i] = 8;
263     ierr = TRUE;
264     break;
265     }
266     }
267     }
268     if (ierr == TRUE)
269     goto L_10;
270     }
271     }
272     /* check four memebered rings */
273     if ( (rings.nring4 > 0) && (angk1.nang4 > 0) )
274     {
275     nRc = is_ring43(ia, ib, ic);
276     if (nRc == TRUE)
277     {
278     /* search four membered ring parameters */
279     for (j = 0; j < angk1.nang4; j++)
280     {
281     if ( (strcmp(angk1.ktype4[j],pt) == 0) || (strcmp(angk1.ktype4[j],pt1) == 0) ||
282     (strcmp(angk1.ktype4[j],pt2) == 0) || (strcmp(angk1.ktype4[j],pt3) == 0) ||
283     (strcmp(angk1.ktype4[j],pt4) == 0) || (strcmp(angk1.ktype4[j],pt5) == 0))
284     {
285     if ( angk1.ang4[j][1] == 0.0 && angk1.ang4[j][2] == 0)
286     {
287     angles.acon[i] = angk1.con4[j];
288     angles.anat[i] = angk1.ang4[j][0];
289     angles.index[i] = 4;
290     ierr = TRUE;
291     break;
292     } else
293     {
294     nh = 0;
295     for (k=0; k < MAXIAT; k++)
296     {
297     ie = atom[ib].iat[k];
298     if ( (atom[ie].type == 5 || atom[ie].type == 36) && (ie != ia) && (ie != ic))
299     nh++;
300     }
301     angles.acon[i] = angk1.con4[j];
302     angles.anat[i] = angk1.ang4[j][nh];
303     angles.index[i] = 4;
304     ierr = TRUE;
305     break;
306     }
307     }
308     }
309     if (ierr == TRUE)
310     goto L_10;
311     else
312     {
313     Missing_constants = TRUE;
314     }
315     }
316     }
317     /* check five membered rings */
318     if ( (rings.nring5 > 0) && (angk1.nang5 > 0) )
319     {
320     nRc = is_ring53(ia,ib, ic);
321     if (nRc == TRUE)
322     {
323     /* search five memebered ring parameters */
324     for (j = 0; j < angk1.nang5; j++)
325     {
326     if ( (strcmp(angk1.ktype5[j],pt) == 0) || (strcmp(angk1.ktype5[j],pt1) == 0) ||
327     (strcmp(angk1.ktype5[j],pt2) == 0) || (strcmp(angk1.ktype5[j],pt3) == 0) ||
328     (strcmp(angk1.ktype5[j],pt4) == 0) || (strcmp(angk1.ktype5[j],pt5) == 0) )
329     {
330     if ( angk1.ang5[j][1] == 0.0 && angk1.ang5[j][2] == 0)
331     {
332     angles.acon[i] = angk1.con5[j];
333     angles.anat[i] = angk1.ang5[j][0];
334     angles.index[i] = 0;
335     ierr = TRUE;
336     break;
337     } else
338     {
339     nh = 0;
340     for (k=0; k < MAXIAT; k++)
341     {
342     ie = atom[ib].iat[k];
343     if ( (atom[ie].type == 5 || atom[ie].type == 36) && (ie != ia) && (ie != ic))
344     nh++;
345     }
346     angles.acon[i] = angk1.con5[j];
347     angles.anat[i] = angk1.ang5[j][nh];
348     angles.index[i] = 0;
349     ierr = TRUE;
350     break;
351     }
352     }
353     }
354     if (ierr == TRUE)
355     goto L_10;
356     // don't fail on missing 5 ring parameters
357     }
358     }
359     /* check delocalized angles in MMFF94 */
360     if (field.type == MMFF94 && angk1.ndel > 0)
361     {
362     if ( bonds_ff.index[nb1] == 1 || bonds_ff.index[nb2] == 1 )
363     {
364     nbo1 = bonds_ff.index[nb1] + bonds_ff.index[nb2];
365     for(j=0; j < angk1.ndel; j++)
366     {
367     if (strcmp(angk1.kdel[j],pt) == 0)
368     {
369     angles.acon[i] = angk1.condel[j];
370     angles.anat[i] = angk1.angdel[j][0];
371     if (nbo1 == 1)
372     angles.index[i] = 1;
373     else if (nbo1 == 2)
374     angles.index[i] = 2;
375     ierr = TRUE;
376     break;
377     }
378     }
379     if (ierr == TRUE)
380     goto L_10;
381     }
382     }
383     // look for fourier angles
384     if (angf.nfang > 0)
385     {
386     for (j=0; j < angf.nfang; j++)
387     {
388     if (strcmp(angf.kftype[j],pt) == 0)
389     {
390     angles.fcon[i] = angf.fcon[j];
391     angles.c0[i] = angf.fc0[j];
392     angles.c1[i] = angf.fc1[j];
393     angles.c2[i] = angf.fc2[j];
394     angles.c3[i] = angf.fc3[j];
395     angles.c4[i] = angf.fc4[j];
396     angles.c5[i] = angf.fc5[j];
397     angles.c6[i] = angf.fc6[j];
398     angles.angtype[i] = FOURIER;
399     ierr = TRUE;
400     break;
401     }
402     }
403     if (ierr == TRUE)
404     goto L_10;
405     }
406    
407     /* check regular parameters for specific angles*/
408     for (j = 0; j < angk1.nang; j++)
409     {
410     if ( (strcmp(angk1.ktype[j],pt) == 0) )
411     {
412     if ( angk1.ang[j][1] == 0.0 && angk1.ang[j][2] == 0)
413     {
414     angles.acon[i] = angk1.con[j];
415     angles.anat[i] = angk1.ang[j][0];
416     angles.index[i] = 0;
417     ierr = TRUE;
418     break;
419     } else
420     {
421     nh = 0;
422     for (k=0; k < MAXIAT; k++)
423     {
424     ie = atom[ib].iat[k];
425     if ( (atom[ie].atomnum == 1) && (ie != ia) && (ie != ic))
426     nh++;
427     }
428     angles.acon[i] = angk1.con[j];
429     angles.anat[i] = angk1.ang[j][nh];
430     angles.index[i] = 0;
431     ierr = TRUE;
432     break;
433     }
434     }
435     }
436     if (ierr == TRUE)
437     goto L_10;
438     // did not find any specific look for generalized
439     for (j = 0; j < angk1.nang; j++)
440     {
441     if ( (strcmp(angk1.ktype[j],pt) == 0) || (strcmp(angk1.ktype[j],pt1) == 0) ||
442     (strcmp(angk1.ktype[j],pt2) == 0) || (strcmp(angk1.ktype[j],pt3) == 0) ||
443     (strcmp(angk1.ktype[j],pt4) == 0) || (strcmp(angk1.ktype[j],pt5) == 0) ||
444     (strcmp(angk1.ktype[j],pt6) == 0))
445     {
446     if ( angk1.ang[j][1] == 0.0 && angk1.ang[j][2] == 0)
447     {
448     angles.acon[i] = angk1.con[j];
449     angles.anat[i] = angk1.ang[j][0];
450     angles.index[i] = 0;
451     ierr = TRUE;
452     break;
453     } else
454     {
455     nh = 0;
456     for (k=0; k < MAXIAT; k++)
457     {
458     ie = atom[ib].iat[k];
459     if ( (atom[ie].atomnum == 1) && (ie != ia) && (ie != ic))
460     nh++;
461     }
462     angles.acon[i] = angk1.con[j];
463     angles.anat[i] = angk1.ang[j][nh];
464     angles.index[i] = 0;
465     ierr = TRUE;
466     break;
467     }
468     }
469     }
470     if (ierr == TRUE)
471     goto L_10;
472     /* check MMFF for class1 parameters */
473     if (field.type == MMFF94)
474     {
475     for (j = 0; j < angk1.nang; j++)
476     {
477     if ( (strcmp(angk1.ktype[j],pt7) == 0) )
478     {
479     if ( angk1.ang[j][1] == 0.0 && angk1.ang[j][2] == 0)
480     {
481     angles.acon[i] = angk1.con[j];
482     angles.anat[i] = angk1.ang[j][0];
483     angles.index[i] = 0;
484     ierr = TRUE;
485     break;
486     } else
487     {
488     nh = 0;
489     for (k=0; k < MAXIAT; k++)
490     {
491     ie = atom[ib].iat[k];
492     if ( (atom[ie].atomnum == 1) && (ie != ia) && (ie != ic))
493     nh++;
494     }
495     angles.acon[i] = angk1.con[j];
496     angles.anat[i] = angk1.ang[j][nh];
497     angles.index[i] = 0;
498     ierr = TRUE;
499     break;
500     }
501     }
502     }
503     if (ierr == TRUE)
504     goto L_10;
505     }
506     if (ierr != TRUE)
507     {
508     bkan = get_general_angle(ibb);
509     angles.acon[i] = 0.80;
510     angles.anat[i] = bkan;
511     angles.index[i] = 0;
512     // ==== done with angle lookup =============
513     L_10:
514     if (field.type == MMX && ierr != TRUE)
515     {
516     if (ibb >= 300)
517     {
518     nligand = 0;
519     jji = 0;
520     for ( k=0; k < MAXIAT; k++)
521     {
522     if (atom[ib].iat[k] != 0 && atom[ib].bo[k] != 9)
523     jji++;
524     if (atom[ib].bo[k] == 9)
525     nligand++;
526     }
527     if (jji == 2 && nligand == 0)
528     {
529     angles.acon[i] = .1;
530     angles.anat[i] = 180.;
531     }
532     if ((jji == 3 && nligand == 0) || (jji == 2 && nligand == 1) )
533     {
534     angles.acon[i] = .1;
535     angles.anat[i] = 120.;
536     }
537     if ((jji == 4 && nligand == 0) || (jji == 3 && nligand == 1) || (jji == 2 && nligand == 2))
538     {
539     if (jji == 4 && (atom[ib].type >= 300 && (atom[ib].flags & (1L << SATMET_MASK) &&
540     atom[ib].flags & (1L << GT18e_MASK) && atom[ib].flags & (1L << SQPLAN_MASK) &&
541     !(atom[ib].flags & (1L << LOWSPIN_MASK)))) )
542     {
543     angles.acon[i] = .0001;
544     angles.anat[i] = 90.;
545     } else
546     {
547     angles.acon[i] = .1;
548     angles.anat[i] = 109.5;
549     }
550     }
551    
552     }
553     }
554     }
555     }
556     }
557     }
558     // ===================================
559     // generalized angles for MMFF with no hydrogen minimization
560     float get_general_angle(int ia)
561     {
562     float angs[80] = {
563     109.0, 120.0, 120.0, 180.0, 0.000, 109.0, 0.000, 109.0, 120.0, 118.0,
564     0.000, 0.000, 0.000, 0.000, 109.0, 109.0, 120.0, 109.0, 109.0, 90.00,
565     0.000, 60.00, 0.000, 0.000, 109.0, 109.0, 0.000, 0.000, 0.000, 90.00,
566     0.000, 0.000, 0.000, 109.0, 0.000, 0.000, 120.0, 120.0, 120.0, 120.0,
567     120.0, 180.0, 120.0, 120.0, 120.0, 120.0, 0.000, 0.000, 109.0, 0.000,
568     109.0, 0.000, 180.0, 120.0, 120.0, 120.0, 120.0, 120.0, 120.0, 180.0,
569     180.0, 120.0, 120.0, 120.0, 120.0, 120.0, 120.0, 109.0, 120.0, 109.0,
570     0.000, 109.0, 109.0, 120.0, 120.0, 109.0, 0.000, 109.0, 109.0, 109.0 };
571    
572     if (ia < 80)
573     return (angs[ia-1]);
574     else
575     return 109.0;
576     }
577     /* --------------------------------- */
578     long kang(int i,int j,int k)
579     {
580     long int kang_v;
581     static long i10000 = 10000;
582     static long i100 = 100;
583    
584     kang_v = i10000*i + i100*j + k;
585     return( kang_v );
586     }
587    
588     /* -------------------- */
589     long ipack(int i,int j,int k,int l)
590     {
591     long int ipack_v;
592     static long itbig = 200000000;
593     static long itmid = 2000000;
594     static long ithou = 1000;
595    
596     ipack_v = itbig*i + itmid*j + ithou*k + l;
597     return( ipack_v );
598     }