ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/kangle.c
Revision: 103
Committed: Thu Feb 19 01:37:38 2009 UTC (12 years, 4 months ago) by gilbertke
File size: 20642 byte(s)
Log Message:
major rewrite - removing global data, adding electrostatics tag to read_sdf
Line File contents
1 #define EXTERN extern
2 #include "pcwin.h"
3 #include "pcmod.h"
4
5 #include "angles.h"
6 #include "bonds_ff.h"
7 #include "atom_k.h"
8
9 int is_ring43(int, int, int);
10 int is_ring53(int, int, int);
11 int is_ring63(int, int, int);
12 int is_delocalbond(int, int);
13 int isbond(int,int);
14 int find_bond(int,int);
15 long ipack(int ,int ,int ,int);
16 long kang(int,int,int);
17 void numeral(int,char *,int);
18 float get_general_angle(int);
19 int have_ring3(void);
20 int have_ring4(void);
21 int have_ring5(void);
22 int have_ring6(void);
23 int get_field(void);
24
25 EXTERN struct t_angk1 {
26 int use_ang3, use_ang4, use_ang5;
27 int nang, nang3, nang4, nang5;
28 int ndel, ndel3, ndel4;
29 char ktype[MAXANGCONST][10],ktype3[MAXANG3CONST][10],ktype4[MAXANG4CONST][10],ktype5[MAXANG5CONST][10];
30 char kdel[MAXANGDEL][10],kdel3[MAXANG3DEL][10],kdel4[MAXANG4DEL][10];
31 float con[MAXANGCONST], ang[MAXANGCONST][3];
32 float con3[MAXANG3CONST], ang3[MAXANG3CONST][3];
33 float con4[MAXANG4CONST], ang4[MAXANG4CONST][3];
34 float con5[MAXANG5CONST], ang5[MAXANG5CONST][3];
35 float condel[MAXANGDEL], angdel[MAXANGDEL][3];
36 float condel3[MAXANG3DEL], angdel3[MAXANG3DEL][3];
37 float condel4[MAXANG4DEL], angdel4[MAXANG4DEL][3];
38 } angk1;
39
40 EXTERN int Missing_constants;
41
42 void kangle()
43 {
44 int i, j, k, izero, ierr, nh, nRc,field;
45 int ia, ib, ic, ie;
46 int iaa, ibb, icc;
47 int cl_a, cl_b, cl_c;
48 int nbo1, nb1, nb2;
49 char pa[4],pb[4],pc[4],pt[10], pt1[10], pt2[10],pt3[10], pt4[10], pt5[10], pt6[10], pt7[10];
50 char zero[4];
51 float bkan;
52
53 if( angles.nang != 0 )
54 {
55 field = get_field();
56 for( i = 0; i < angles.nang; i++ )
57 {
58 ia = angles.i13[i][0];
59 ib = angles.i13[i][1];
60 ic = angles.i13[i][2];
61 iaa = atom.type[ia];
62 ibb = atom.type[ib];
63 icc = atom.type[ic];
64 cl_a = atom.tclass[ia];
65 cl_b = atom.tclass[ib];
66 cl_c = atom.tclass[ic];
67
68 /* get bond index each bond for MMFF94 */
69 nb1 = find_bond(ia,ib);
70 nb2 = find_bond(ib,ic);
71
72 /* special metal cases */
73 if( iaa >= 300)
74 cl_a = 300;
75
76 if( ibb >= 300)
77 cl_b = 300;
78
79 if( icc >= 300)
80 cl_c = 300;
81
82 /* special MMX cases */
83 if ( field == MMX)
84 {
85 if( iaa == 40 )
86 iaa = 2;
87
88 if( ibb == 40 )
89 ibb = 2;
90
91 if( icc == 40 )
92 icc = 2;
93
94 }
95 numeral(iaa,pa,3);
96 numeral(ibb,pb,3);
97 numeral(icc,pc,3);
98 if( iaa <= icc )
99 {
100 strcpy(pt,pa);
101 strcat(pt,pb);
102 strcat(pt,pc);
103 }
104 if( iaa > icc )
105 {
106 strcpy(pt,pc);
107 strcat(pt,pb);
108 strcat(pt,pa);
109 }
110 izero = 0;
111 strcpy(zero," 0");
112 strcpy(pt1,pa); strcat(pt1,pb); strcat(pt1,zero);
113 strcpy(pt2,pc); strcat(pt2,pb); strcat(pt2,zero);
114 strcpy(pt3,zero); strcat(pt3,pb); strcat(pt3,pc);
115 strcpy(pt4,zero); strcat(pt4,pb); strcat(pt4,pa);
116 strcpy(pt5,zero); strcat(pt5,pb); strcat(pt5,zero);
117
118 numeral(cl_a,pa,3);
119 numeral(cl_b,pb,3);
120 numeral(cl_c,pc,3);
121 if (cl_a < cl_c)
122 {
123 strcpy(pt6,pa); strcat(pt6,pb); strcat(pt6,pc);
124 }else
125 {
126 strcpy(pt6,pc); strcat(pt6,pb); strcat(pt6,pa);
127 }
128
129 if (field == MMFF94)
130 {
131 cl_a = atom_k.tclass1[atom.type[ia]];
132 cl_b = atom_k.tclass1[atom.type[ib]];
133 cl_c = atom_k.tclass1[atom.type[ic]];
134 numeral(cl_a,pa,3);
135 numeral(cl_b,pb,3);
136 numeral(cl_c,pc,3);
137 if (cl_a < cl_c)
138 {
139 strcpy(pt7,pa); strcat(pt7,pb); strcat(pt7,pc);
140 }else
141 {
142 strcpy(pt7,pc); strcat(pt7,pb); strcat(pt7,pa);
143 }
144 }
145 angles.acon[i] = 0.0;
146 angles.anat[i] = 0.0;
147 angles.angtype[i] = HARMONIC;
148
149 ierr = FALSE;
150 /* check TS constants */
151 /* check delocalized angles in MMFF94 cyclopropanes */
152 if (field == MMFF94 && angk1.ndel3 > 0)
153 {
154 if (isbond(ia,ic) )
155 {
156 nbo1 = bonds_ff.index[nb1] + bonds_ff.index[nb2];
157 if (nbo1 == 1 || nbo1 == 2)
158 {
159 for(j=0; j < angk1.ndel3; j++)
160 {
161 if (strcmp(angk1.kdel3[j],pt) == 0)
162 {
163 angles.acon[i] = angk1.condel3[j];
164 angles.anat[i] = angk1.angdel3[j][0];
165 if (nbo1 == 1)
166 angles.index[i] = 5;
167 else if (nbo1 == 2)
168 angles.index[i] = 6;
169 ierr = TRUE;
170 break;
171 }
172 }
173 }
174
175 if (ierr == TRUE)
176 goto L_10;
177 }
178 }
179 /* check three membered rings */
180 if ( (have_ring3()) && (angk1.nang3 > 0) )
181 {
182 if (isbond(ia,ic))
183 {
184 /* search three membered ring parameters */
185 for (j = 0; j < angk1.nang3; j++)
186 {
187 if ( (strcmp(angk1.ktype3[j],pt) == 0) || (strcmp(angk1.ktype3[j],pt1) == 0) ||
188 (strcmp(angk1.ktype3[j],pt2) == 0) || (strcmp(angk1.ktype3[j],pt3) == 0) ||
189 (strcmp(angk1.ktype3[j],pt4) == 0) || (strcmp(angk1.ktype3[j],pt5) == 0))
190 {
191 if ( angk1.ang3[j][1] == 0.0 && angk1.ang3[j][2] == 0)
192 {
193 angles.acon[i] = angk1.con3[j];
194 angles.anat[i] = angk1.ang3[j][0];
195 angles.index[i] = 3;
196 ierr = TRUE;
197 break;
198 } else
199 {
200 nh = 0;
201 for (k=0; k < MAXIAT; k++)
202 {
203 ie = atom.iat[ib][k];
204 if ( (atom.type[ie] == 5 || atom.type[ie] == 36) && (ie != ia) && (ie != ic))
205 nh++;
206 }
207 angles.acon[i] = angk1.con3[j];
208 angles.anat[i] = angk1.ang3[j][nh];
209 angles.index[i] = 3;
210 ierr = TRUE;
211 break;
212 }
213 }
214 }
215 if (ierr == TRUE)
216 goto L_10;
217 else
218 {
219 Missing_constants = TRUE;
220 }
221 }
222 }
223 /* check delocalized angles in MMFF94 cyclobutanes */
224 if (field == MMFF94 && angk1.ndel4 > 0)
225 {
226 nRc = is_ring43(ia, ib, ic);
227 if (nRc == TRUE)
228 {
229 nbo1 = bonds_ff.index[nb1] + bonds_ff.index[nb2];
230 if (nbo1 == 1 || nbo1 == 2)
231 {
232 for(j=0; j < angk1.ndel4; j++)
233 {
234 if (strcmp(angk1.kdel4[j],pt) == 0)
235 {
236 angles.acon[i] = angk1.condel4[j];
237 angles.anat[i] = angk1.angdel4[j][0];
238 if (nbo1 == 1)
239 angles.index[i] = 7;
240 else if (nbo1 == 2)
241 angles.index[i] = 8;
242 ierr = TRUE;
243 break;
244 }
245 }
246 }
247 if (ierr == TRUE)
248 goto L_10;
249 }
250 }
251 /* check four memebered rings */
252 if ( (have_ring4()) && (angk1.nang4 > 0) )
253 {
254 nRc = is_ring43(ia, ib, ic);
255 if (nRc == TRUE)
256 {
257 /* search four membered ring parameters */
258 for (j = 0; j < angk1.nang4; j++)
259 {
260 if ( (strcmp(angk1.ktype4[j],pt) == 0) || (strcmp(angk1.ktype4[j],pt1) == 0) ||
261 (strcmp(angk1.ktype4[j],pt2) == 0) || (strcmp(angk1.ktype4[j],pt3) == 0) ||
262 (strcmp(angk1.ktype4[j],pt4) == 0) || (strcmp(angk1.ktype4[j],pt5) == 0))
263 {
264 if ( angk1.ang4[j][1] == 0.0 && angk1.ang4[j][2] == 0)
265 {
266 angles.acon[i] = angk1.con4[j];
267 angles.anat[i] = angk1.ang4[j][0];
268 angles.index[i] = 4;
269 ierr = TRUE;
270 break;
271 } else
272 {
273 nh = 0;
274 for (k=0; k < MAXIAT; k++)
275 {
276 ie = atom.iat[ib][k];
277 if ( (atom.type[ie] == 5 || atom.type[ie] == 36) && (ie != ia) && (ie != ic))
278 nh++;
279 }
280 angles.acon[i] = angk1.con4[j];
281 angles.anat[i] = angk1.ang4[j][nh];
282 angles.index[i] = 4;
283 ierr = TRUE;
284 break;
285 }
286 }
287 }
288 if (ierr == TRUE)
289 goto L_10;
290 else
291 {
292 Missing_constants = TRUE;
293 }
294 }
295 }
296 /* check five membered rings */
297 if ( (have_ring5()) && (angk1.nang5 > 0) )
298 {
299 nRc = is_ring53(ia,ib, ic);
300 if (nRc == TRUE)
301 {
302 /* search five memebered ring parameters */
303 for (j = 0; j < angk1.nang5; j++)
304 {
305 if ( (strcmp(angk1.ktype5[j],pt) == 0) || (strcmp(angk1.ktype5[j],pt1) == 0) ||
306 (strcmp(angk1.ktype5[j],pt2) == 0) || (strcmp(angk1.ktype5[j],pt3) == 0) ||
307 (strcmp(angk1.ktype5[j],pt4) == 0) || (strcmp(angk1.ktype5[j],pt5) == 0) )
308 {
309 if ( angk1.ang5[j][1] == 0.0 && angk1.ang5[j][2] == 0)
310 {
311 angles.acon[i] = angk1.con5[j];
312 angles.anat[i] = angk1.ang5[j][0];
313 angles.index[i] = 0;
314 ierr = TRUE;
315 break;
316 } else
317 {
318 nh = 0;
319 for (k=0; k < MAXIAT; k++)
320 {
321 ie = atom.iat[ib][k];
322 if ( (atom.type[ie] == 5 || atom.type[ie] == 36) && (ie != ia) && (ie != ic))
323 nh++;
324 }
325 angles.acon[i] = angk1.con5[j];
326 angles.anat[i] = angk1.ang5[j][nh];
327 angles.index[i] = 0;
328 ierr = TRUE;
329 break;
330 }
331 }
332 }
333 if (ierr == TRUE)
334 goto L_10;
335 // don't fail on missing 5 ring parameters
336 }
337 }
338 /* check delocalized angles in MMFF94 */
339 if (field == MMFF94 && angk1.ndel > 0)
340 {
341 if ( bonds_ff.index[nb1] == 1 || bonds_ff.index[nb2] == 1 )
342 {
343 nbo1 = bonds_ff.index[nb1] + bonds_ff.index[nb2];
344 for(j=0; j < angk1.ndel; j++)
345 {
346 if (strcmp(angk1.kdel[j],pt) == 0)
347 {
348 angles.acon[i] = angk1.condel[j];
349 angles.anat[i] = angk1.angdel[j][0];
350 if (nbo1 == 1)
351 angles.index[i] = 1;
352 else if (nbo1 == 2)
353 angles.index[i] = 2;
354 ierr = TRUE;
355 break;
356 }
357 }
358 if (ierr == TRUE)
359 goto L_10;
360 }
361 }
362 /* check regular parameters for specific angles*/
363 for (j = 0; j < angk1.nang; j++)
364 {
365 if ( (strcmp(angk1.ktype[j],pt) == 0) )
366 {
367 if ( angk1.ang[j][1] == 0.0 && angk1.ang[j][2] == 0)
368 {
369 angles.acon[i] = angk1.con[j];
370 angles.anat[i] = angk1.ang[j][0];
371 angles.index[i] = 0;
372 ierr = TRUE;
373 break;
374 } else
375 {
376 nh = 0;
377 for (k=0; k < MAXIAT; k++)
378 {
379 ie = atom.iat[ib][k];
380 if ( (atom.atomnum[ie] == 1) && (ie != ia) && (ie != ic))
381 nh++;
382 }
383 angles.acon[i] = angk1.con[j];
384 angles.anat[i] = angk1.ang[j][nh];
385 angles.index[i] = 0;
386 ierr = TRUE;
387 break;
388 }
389 }
390 }
391 if (ierr == TRUE)
392 goto L_10;
393 // did not find any specific look for generalized
394 for (j = 0; j < angk1.nang; j++)
395 {
396 if ( (strcmp(angk1.ktype[j],pt) == 0) || (strcmp(angk1.ktype[j],pt1) == 0) ||
397 (strcmp(angk1.ktype[j],pt2) == 0) || (strcmp(angk1.ktype[j],pt3) == 0) ||
398 (strcmp(angk1.ktype[j],pt4) == 0) || (strcmp(angk1.ktype[j],pt5) == 0) ||
399 (strcmp(angk1.ktype[j],pt6) == 0))
400 {
401 if ( angk1.ang[j][1] == 0.0 && angk1.ang[j][2] == 0)
402 {
403 angles.acon[i] = angk1.con[j];
404 angles.anat[i] = angk1.ang[j][0];
405 angles.index[i] = 0;
406 ierr = TRUE;
407 break;
408 } else
409 {
410 nh = 0;
411 for (k=0; k < MAXIAT; k++)
412 {
413 ie = atom.iat[ib][k];
414 if ( (atom.atomnum[ie] == 1) && (ie != ia) && (ie != ic))
415 nh++;
416 }
417 angles.acon[i] = angk1.con[j];
418 angles.anat[i] = angk1.ang[j][nh];
419 angles.index[i] = 0;
420 ierr = TRUE;
421 break;
422 }
423 }
424 }
425 if (ierr == TRUE)
426 goto L_10;
427 /* check MMFF for class1 parameters */
428 if (field == MMFF94)
429 {
430 for (j = 0; j < angk1.nang; j++)
431 {
432 if ( (strcmp(angk1.ktype[j],pt7) == 0) )
433 {
434 if ( angk1.ang[j][1] == 0.0 && angk1.ang[j][2] == 0)
435 {
436 angles.acon[i] = angk1.con[j];
437 angles.anat[i] = angk1.ang[j][0];
438 angles.index[i] = 0;
439 ierr = TRUE;
440 break;
441 } else
442 {
443 nh = 0;
444 for (k=0; k < MAXIAT; k++)
445 {
446 ie = atom.iat[ib][k];
447 if ( (atom.atomnum[ie] == 1) && (ie != ia) && (ie != ic))
448 nh++;
449 }
450 angles.acon[i] = angk1.con[j];
451 angles.anat[i] = angk1.ang[j][nh];
452 angles.index[i] = 0;
453 ierr = TRUE;
454 break;
455 }
456 }
457 }
458 if (ierr == TRUE)
459 goto L_10;
460 }
461 if (ierr != TRUE)
462 {
463 bkan = get_general_angle(ibb);
464 angles.acon[i] = 0.80;
465 angles.anat[i] = bkan;
466 angles.index[i] = 0;
467 }
468 // ==== done with angle lookup =============
469 L_10:
470 continue;
471 }
472 }
473 }
474 // ===================================
475 // generalized angles for MMFF with no hydrogen minimization
476 float get_general_angle(int ia)
477 {
478 float angs[80] = {
479 109.0, 120.0, 120.0, 180.0, 0.000, 109.0, 0.000, 109.0, 120.0, 118.0,
480 0.000, 0.000, 0.000, 0.000, 109.0, 109.0, 120.0, 109.0, 109.0, 90.00,
481 0.000, 60.00, 0.000, 0.000, 109.0, 109.0, 0.000, 0.000, 0.000, 90.00,
482 0.000, 0.000, 0.000, 109.0, 0.000, 0.000, 120.0, 120.0, 120.0, 120.0,
483 120.0, 180.0, 120.0, 120.0, 120.0, 120.0, 0.000, 0.000, 109.0, 0.000,
484 109.0, 0.000, 180.0, 120.0, 120.0, 120.0, 120.0, 120.0, 120.0, 180.0,
485 180.0, 120.0, 120.0, 120.0, 120.0, 120.0, 120.0, 109.0, 120.0, 109.0,
486 0.000, 109.0, 109.0, 120.0, 120.0, 109.0, 0.000, 109.0, 109.0, 109.0 };
487
488 if (ia < 80)
489 return (angs[ia-1]);
490 else
491 return 109.0;
492 }
493 /* --------------------------------- */
494 long kang(int i,int j,int k)
495 {
496 long int kang_v;
497 static long i10000 = 10000;
498 static long i100 = 100;
499
500 kang_v = i10000*i + i100*j + k;
501 return( kang_v );
502 }
503
504 /* -------------------- */
505 long ipack(int i,int j,int k,int l)
506 {
507 long int ipack_v;
508 static long itbig = 200000000;
509 static long itmid = 2000000;
510 static long ithou = 1000;
511
512 ipack_v = itbig*i + itmid*j + ithou*k + l;
513 return( ipack_v );
514 }