ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/kbond.c
Revision: 64
Committed: Thu Dec 4 02:03:23 2008 UTC (13 years, 1 month ago) by gilbertke
File size: 19243 byte(s)
Log Message:
fixed bug in routine to recognize delocalized bonds in kbond.c
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5 #include "bonds_ff.h"
6 #include "rings.h"
7 #include "field.h"
8
9 FILE * fopen_path ( char * , char * , char * ) ;
10 void angle(int ,int ,int ,float *);
11 int is_ring31(int);
12 int is_ring41(int);
13 int is_ring51(int);
14 int is_ring61(int);
15 int is_ring42(int, int);
16 int is_ring52(int, int);
17 int is_ring62(int,int);
18 void pcmfout(int);
19 double arcos(float );
20 double angs(float,float,float,float,float ,float ,float ,float);
21 double powi(double, int);
22 void numeral(int, char *, int);
23 int is_delocalbond(int ia, int ib);
24 void message_alert(char *, char *);
25 double SIGN(double ,double );
26 int is_bond(int,int);
27 double deltaks(double,double);
28 float get_bond(int, int);
29
30
31 EXTERN struct t_bondk1 {
32 int use_ring3, use_ring4, use_ring5;
33 int nbnd, nbnd3, nbnd4, nbnd5, ndeloc;
34 char kb[MAXBONDCONST][7], kb3[MAXBOND3CONST][7],kb4[MAXBOND4CONST][7],kb5[MAXBOND5CONST][7];
35 char kbdel[MAXBONDDELOC][7];
36 float s[MAXBONDCONST], t[MAXBONDCONST];
37 float s3[MAXBOND3CONST], t3[MAXBOND3CONST];
38 float s4[MAXBOND4CONST], t4[MAXBOND4CONST];
39 float s5[MAXBOND5CONST], t5[MAXBOND5CONST];
40 float sdel[MAXBONDDELOC], tdel[MAXBONDDELOC];
41 } bondk1;
42
43
44 EXTERN struct t_bondpk {
45 int npibond,npibond44i,npibond45i,npibond46i,npibond55i,npibond56i,npibond66i;
46 int npibond44o,npibond45o,npibond46o,npibond55o,npibond56o,npibond66o;
47 int npibond40,npibond50,npibond60;
48 int use_pibond,use_pibond44i,use_pibond45i,use_pibond46i,use_pibond55i,use_pibond56i,use_pibond66i;
49 int use_pibond44o,use_pibond45o,use_pibond46o,use_pibond55o,use_pibond56o,use_pibond66o;
50 int use_pibond40,use_pibond50,use_pibond60;
51 char kb[MAXPIBONDCONST][7],kb44i[20][7],kb45i[20][7],kb46i[20][7],kb55i[20][7],kb56i[20][7],kb66i[20][7];
52 char kb44o[20][7],kb45o[20][7],kb46o[20][7],kb55o[20][7],kb56o[20][7],kb66o[20][7];
53 char kb40[20][7],kb50[20][7],kb60[20][7];
54 float bk[MAXPIBONDCONST], bl[MAXPIBONDCONST], bmom[MAXPIBONDCONST],sslop[MAXPIBONDCONST], tslop[MAXPIBONDCONST], tslop2[MAXPIBONDCONST];
55 float bk44i[20], bl44i[20], bmom44i[20],sslop44i[20], tslop44i[20], tslop244i[20];
56 float bk45i[20], bl45i[20], bmom45i[20],sslop45i[20], tslop45i[20], tslop245i[20];
57 float bk46i[20], bl46i[20], bmom46i[20],sslop46i[20], tslop46i[20], tslop246i[20];
58 float bk55i[20], bl55i[20], bmom55i[20],sslop55i[20], tslop55i[20], tslop255i[20];
59 float bk56i[20], bl56i[20], bmom56i[20],sslop56i[20], tslop56i[20], tslop256i[20];
60 float bk66i[20], bl66i[20], bmom66i[20],sslop66i[20], tslop66i[20], tslop266i[20];
61 float bk44o[20], bl44o[20], bmom44o[20],sslop44o[20], tslop44o[20], tslop244o[20];
62 float bk45o[20], bl45o[20], bmom45o[20],sslop45o[20], tslop45o[20], tslop245o[20];
63 float bk46o[20], bl46o[20], bmom46o[20],sslop46o[20], tslop46o[20], tslop246o[20];
64 float bk55o[20], bl55o[20], bmom55o[20],sslop55o[20], tslop55o[20], tslop255o[20];
65 float bk56o[20], bl56o[20], bmom56o[20],sslop56o[20], tslop56o[20], tslop256o[20];
66 float bk66o[20], bl66o[20], bmom66o[20],sslop66o[20], tslop66o[20], tslop266o[20];
67 float bk40[20], bl40[20], bmom40[20],sslop40[20], tslop40[20], tslop240[20];
68 float bk50[20], bl50[20], bmom50[20],sslop50[20], tslop50[20], tslop250[20];
69 float bk60[20], bl60[20], bmom60[20],sslop60[20], tslop60[20], tslop260[20];
70 } bondpk;
71
72 EXTERN struct t_electroneg {
73 int nelecti, nelectj;
74 int itype[200],ibond[200],iattach[200];
75 int jtype[200],jbond[200],jattach[200];
76 float icorr[200],jcorr[200];
77 } electroneg;
78
79 EXTERN struct t_files {
80 int nfiles, append, batch, icurrent;
81 int istrselect;
82 } files;
83
84 EXTERN int Missing_constants;
85 EXTERN FILE *errfile;
86
87 int kbond()
88 {
89 long int pi_mask;
90 int ia4,ia5,ia6, ib4,ib5,ib6, iab4, iab5, iab6;
91 int i, j,k, iit, kit, it, ierr,nbpi, ii, kk, i1, i2;
92 int igeni, igenk, ifound;
93 int ia, ib, iend, kend, iy;
94 int numi, numj, iatt_electro[6], jatt_electro[6];
95 float jatt_const[6],iatt_const[6];
96 float bconst, blength;
97 float sslope, tslope, tslope2;
98 double dbl,dks,bl0;
99 char hatext[80];
100 char pa[4],pb[4], pt[7], pt_gen[7];
101
102 pi_mask = 1L << PI_MASK; /* this mask is for pi atoms - now uses bit 0 */
103 nbpi = 0;
104 if( natom != 0 )
105 {
106
107 for( i = 0; i < bonds_ff.nbnd; i++ )
108 {
109 ia = bonds_ff.i12[i][0];
110 ib = bonds_ff.i12[i][1];
111 iit = atom[bonds_ff.i12[i][0]].type;
112 kit = atom[bonds_ff.i12[i][1]].type;
113
114 igeni = atom[bonds_ff.i12[i][0]].tclass;
115 igenk = atom[bonds_ff.i12[i][1]].tclass;
116
117 // metal general class
118 if (iit > 300)
119 igeni = 300;
120 if (kit > 300)
121 igenk = 300;
122
123 if( field.type == MMX)
124 {
125 if (iit == 40 && ( kit != 2 && kit != 3 && kit != 4 && kit != 40) )
126 iit = 2;
127 if (kit == 40 && ( iit != 2 && iit != 3 && iit != 4 && iit != 40) )
128 kit = 2;
129 }
130
131 numeral(iit, pa, 3);
132 numeral(kit, pb, 3);
133 strcpy(pt,pa);
134 strcat(pt,pb);
135 if( iit <= kit )
136 {
137 strcpy(pt,pa);
138 strcat(pt,pb);
139 }else
140 {
141 strcpy(pt,pb);
142 strcat(pt,pa);
143 }
144
145 numeral(igeni, pa, 3);
146 numeral(igenk, pb, 3);
147 if( igeni <= igenk )
148 {
149 strcpy(pt_gen,pa);
150 strcat(pt_gen,pb);
151 }else
152 {
153 strcpy(pt_gen,pb);
154 strcat(pt_gen,pa);
155 }
156
157 ierr = FALSE;
158 // check 3 membered ring
159 if (rings.nring3 > 0 && bondk1.nbnd3 > 0 && ( is_ring31(ia) && is_ring31(ib) ) )
160 {
161 for (j=0; j < bondk1.nbnd3; j++)
162 {
163 if ( (strcmp(bondk1.kb3[j], pt) == 0) || (strcmp(bondk1.kb3[j],pt_gen) == 0) )
164 {
165 bonds_ff.bk[i] = bondk1.s3[j];
166 bonds_ff.bl[i] = bondk1.t3[j];
167 bonds_ff.index[i] = 0;
168 ierr = TRUE;
169 break;
170 }
171 }
172 if (ierr != TRUE)
173 {
174 Missing_constants = TRUE;
175 sprintf(hatext,"Bond constants missing for cyclopropane bond %d - %d of type %d-%d\n",bonds_ff.i12[i][0],bonds_ff.i12[i][1],
176 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
177 fprintf(pcmlogfile,"%s\n",hatext);
178 //printf("%s\n",hatext);
179 }
180 }
181 if (ierr == TRUE)
182 goto L_10;
183 // check 4 membered ring
184 if (rings.nring4 > 0 && bondk1.nbnd4 > 0 && ( is_ring41(ia) && is_ring41(ib) ) )
185 {
186 for (j=0; j < bondk1.nbnd4; j++)
187 {
188 if ( (strcmp(bondk1.kb4[j],pt) == 0) || (strcmp(bondk1.kb4[j],pt_gen) == 0) )
189 {
190 bonds_ff.bk[i] = bondk1.s4[j];
191 bonds_ff.bl[i] = bondk1.t4[j];
192 bonds_ff.index[i] = 0;
193 ierr = TRUE;
194 break;
195 }
196 }
197 if (ierr != TRUE)
198 {
199 Missing_constants = TRUE;
200 sprintf(hatext,"Bond constants missing for cyclobutane bond %d - %d of type %d-%d\n",bonds_ff.i12[i][0],bonds_ff.i12[i][1],
201 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
202 fprintf(pcmlogfile,"%s\n",hatext);
203 // printf("%s\n",hatext);
204 }
205 }
206 if (ierr == TRUE)
207 goto L_10;
208 // check 5 membered ring
209 if (rings.nring5 > 0 && bondk1.nbnd5 > 0 && (is_ring51(ia) != -1) && (is_ring51(ib) != -1))
210 {
211 for (j=0; j < bondk1.nbnd5; j++)
212 {
213 if (strcmp(bondk1.kb5[j],pt) == 0 )
214 {
215 bonds_ff.bk[i] = bondk1.s5[j];
216 bonds_ff.bl[i] = bondk1.t5[j];
217 bonds_ff.index[i] = 0;
218 ierr = TRUE;
219 break;
220 }
221 }
222 // don't fail on missing param - check for regular param first
223 }
224 if (ierr == TRUE)
225 goto L_10;
226 // delocalized bonds in MMFF94
227 if (field.type == MMFF94 && bondk1.ndeloc > 0)
228 {
229 if ( is_delocalbond(ia, ib) )
230 {
231 for (k = 0; k < MAXIAT; k++)
232 {
233 if (atom[ia].iat[k] == ib)
234 {
235 if (atom[ia].bo[k] == 1) // single bond
236 {
237 for (j=0; j < bondk1.ndeloc; j++)
238 {
239 if (strcmp(bondk1.kbdel[j],pt) == 0)
240 {
241 bonds_ff.bk[i] = bondk1.sdel[j];
242 bonds_ff.bl[i] = bondk1.tdel[j];
243 bonds_ff.index[i] = 1;
244 ierr = TRUE;
245 break;
246 }
247 }
248 }
249 }
250 }
251 }
252 if (ierr == TRUE)
253 goto L_10;
254 }
255 // regular parameters
256 for (j=0; j < bondk1.nbnd; j++)
257 {
258 if ( (strcmp(bondk1.kb[j], pt) == 0) )
259 {
260 bonds_ff.bk[i] = bondk1.s[j];
261 bonds_ff.bl[i] = bondk1.t[j];
262 bonds_ff.index[i] = 0;
263 ierr = TRUE;
264 if ( ((strcmp(pt," 2 2") == 0) || (strcmp(pt," 2 3") == 0))&& field.type == MMX)
265 {
266 for (k=0; k < MAXIAT; k++)
267 {
268 if (atom[ia].iat[k] == ib)
269 {
270 if (atom[ia].bo[k] == 1) // butadiene
271 {
272 bonds_ff.bl[i] = 1.48;
273 break;
274 }
275 }
276 }
277 }
278 break;
279 }
280 }
281 if (ierr == TRUE)
282 goto L_10;
283 //
284 // generalized parameters
285 for (j=0; j < bondk1.nbnd; j++)
286 {
287 if ( (strcmp(bondk1.kb[j], pt_gen) == 0) )
288 {
289 bonds_ff.bk[i] = bondk1.s[j];
290 bonds_ff.bl[i] = bondk1.t[j];
291 bonds_ff.index[i] = 0;
292 ierr = TRUE;
293 break;
294 }
295 }
296 if (ierr == TRUE)
297 goto L_10;
298 L_10:
299 continue;
300
301 }
302 }
303 if (field.type == MM3) // need check for electronegativity correction
304 {
305 for( i = 0; i < bonds_ff.nbnd; i++ )
306 {
307 ia = bonds_ff.i12[i][0];
308 ib = bonds_ff.i12[i][1];
309 iit = atom[bonds_ff.i12[i][0]].tclass;
310 kit = atom[bonds_ff.i12[i][1]].tclass;
311 if (iit > kit)
312 {
313 iend = ib;
314 kend = ia;
315 it = 1000*kit+iit;
316 }else
317 {
318 iend = ia;
319 kend = ib;
320 it = 1000*iit+kit;
321 }
322 numi = 0;
323 for (j=0; j < 6; j++)
324 {
325 iatt_electro[j] = 0;
326 iatt_const[j] = 0.0;
327 jatt_electro[j] = 0;
328 jatt_const[j] = 0.0;
329 }
330 for (j=0; j < electroneg.nelecti; j++)
331 {
332 if (it == electroneg.ibond[j])
333 {
334 for (k=0; k < MAXIAT; k++)
335 {
336 if (atom[iend].iat[k] != kend && atom[atom[iend].iat[k]].tclass == electroneg.iattach[j])
337 {
338 iatt_electro[numi]++;
339 iatt_const[numi] = electroneg.icorr[j];
340 numi++;
341 if (numi > 5) message_alert("error in kbond numi > 5","Error");
342 }
343 }
344 if (iit == kit)
345 {
346 for (k=0; k < MAXIAT; k++)
347 {
348 if (atom[kend].iat[k] != iend && atom[atom[kend].iat[k]].tclass == electroneg.iattach[j])
349 {
350 iatt_electro[numi]++;
351 iatt_const[numi] = electroneg.icorr[j];
352 numi++;
353 if (numi > 5) message_alert("error in kbond numi > 5","Error");
354 }
355 }
356 }
357 }
358 }
359 // loop over ib " " " "
360 numj = 0;
361 jatt_electro[0] = jatt_electro[1] = jatt_electro[2] = 0;
362 jatt_const[0] = jatt_const[1] = jatt_const[2] = 0.0;
363 for (j=0; j < electroneg.nelectj; j++)
364 {
365 if (it == electroneg.jbond[j])
366 {
367 for (k=0; k < MAXIAT; k++)
368 {
369 if (atom[kend].iat[k] != iend && atom[atom[kend].iat[k]].tclass == electroneg.jattach[j])
370 {
371 jatt_electro[numj]++;
372 jatt_const[numj] = electroneg.jcorr[j];
373 numj++;
374 if (numj > 5) message_alert("error in kbond numj > 5","Error");
375 }
376 }
377 }
378 }
379 // correct bond length
380 dbl = 0.0;
381 iy = -1;
382 for (j=0; j < numi; j++)
383 {
384 iy++;
385 dbl += iatt_const[j]*powi(0.62,iy);
386 }
387 for (j=0; j < numj; j++)
388 {
389 iy++;
390 dbl += jatt_const[j]*powi(0.62,iy);
391 }
392 bonds_ff.bl[i] += dbl;
393 //adjust force const
394 if (dbl != 0.0)
395 {
396 if (it == 1005 || it == 5022 || it == 5056)
397 {
398 bl0 = bonds_ff.bl[i] - dbl;
399 dks = deltaks(dbl,bl0);
400 bonds_ff.bk[i] += dks;
401 }
402 }
403
404 }
405 }
406 if (Missing_constants == TRUE)
407 {
408 fprintf(pcmlogfile,"Error assigning constants in: %s\n",Struct_Title);
409 fprintf(pcmlogfile,"%s\n",hatext);
410 return FALSE;
411 }
412
413 return TRUE;
414 }
415 // ==============================================
416 double deltaks(double dbl,double bl0)
417 {
418 double a,b,pi,clite,amu,conv,delta;
419 a = 1.3982;
420 b = -0.0001023;
421 pi = acos(-1.0);
422 clite = 2.9979e10;
423 amu = 1.660531e-24;
424 conv = 4.0*pi*pi*clite*clite*1.008*amu*1.0e-5;
425 delta = conv*dbl*2.0*(bl0-a)/(b*b);
426 return(delta);
427 }
428 /* ------------------------------------------ */
429 int find_bond(int ia, int ib)
430 {
431 int i;
432 int iz;
433
434 iz = -1;
435 for (i=0; i < bonds_ff.nbnd; i++)
436 {
437 if (ia == bonds_ff.i12[i][0] && ib == bonds_ff.i12[i][1])
438 return(i);
439 if (ib == bonds_ff.i12[i][0] && ia == bonds_ff.i12[i][1])
440 return(i);
441 }
442 return(iz);
443 }
444
445 /* ------------------------------------------ */
446 void angle(int k,int i,int imi,float *a13)
447 {
448 float disik, disim, dx1,dy1,dz1, dx2,dy2,dz2;
449
450 /* returns angle between k-i-imi */
451
452 dx1 = atom[i].x - atom[k].x;
453 dy1 = atom[i].y - atom[k].y;
454 dz1 = atom[i].z - atom[k].z;
455 dx2 = atom[i].x - atom[imi].x;
456 dy2 = atom[i].y - atom[imi].y;
457 dz2 = atom[i].z - atom[imi].z;
458
459 disik = sqrt( dx1*dx1 + dy1*dy1 + dz1*dz1);
460 disim = sqrt( dx2*dx2 + dy2*dy2 + dz2*dz2);
461
462 *a13 = angs( dx1, dy1, dz1, dx2, dy2, dz2, disik, disim );
463
464 *a13 *= 57.2958;
465
466 return;
467 }
468 /* ------------------------------ */
469 double angs(float x1,float y1,float z1,float x2,float y2,float z2,float d1,float d2)
470 {
471 double angs_v, cosa;
472
473 /* *** calculates angles for subroutine ebend *** */
474 cosa = (x1*x2 + y1*y2 + z1*z2)/(d1*d2);
475 if( fabs( cosa ) > 1.0 )
476 cosa /= fabs( cosa );
477 angs_v = arcos( cosa );
478 return( angs_v );
479 }
480 /* ----------------------- */
481 double arcos(float x)
482 {
483 double arcos_v, y;
484
485 y = x;
486 if( y > 1.0 )
487 y = 1.0;
488 if( y < -1.0 )
489 y = -1.0;
490 arcos_v = acos( y );
491 return( arcos_v );
492 }
493 /* ------------------------------------------ */
494 int is_delocalbond(int ia, int ib)
495 {
496 // test for delocalized bond
497 // if ia & ib are aromatic and part of 6 membered ring
498 // if not aromatic is bond order 1
499
500 long int aromatic_mask;
501 int i, j;
502 int jdbl, kdbl;
503
504 aromatic_mask = (1L << AROMATIC_MASK);
505
506 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
507 && (is_ring62(ia,ib) != FALSE) )
508 return(FALSE);
509
510 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
511 && (is_ring52(ia,ib) != FALSE) )
512 return(FALSE);
513
514 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask) )
515 return(TRUE);
516 if (atom[ia].type == 37 && atom[ib].type == 37 && (is_ring62(ia,ib) != FALSE))
517 return FALSE;
518 if (atom[ia].type == 37 && atom[ib].type == 37 && (is_ring52(ia,ib) != FALSE))
519 return FALSE;
520
521 for (i=0; i < MAXIAT; i++)
522 {
523 if (atom[ia].iat[i] == ib)
524 {
525 if (atom[ia].bo[i] == 1)
526 {
527 jdbl = FALSE;
528 kdbl = FALSE;
529 for (j=0; j < MAXIAT; j++)
530 {
531 if (atom[ia].bo[j] >= 2 && atom[ia].bo[j] != 9)
532 jdbl = TRUE;
533 if (atom[ib].bo[j] >= 2 && atom[ib].bo[j] != 9)
534 kdbl = TRUE;
535 }
536 if (jdbl == TRUE && kdbl == TRUE)
537 return (TRUE);
538 else if (jdbl == TRUE && (atom[ib].flags & aromatic_mask) )
539 return (TRUE);
540 else if (kdbl == TRUE && (atom[ia].flags & aromatic_mask) )
541 return (TRUE);
542 else
543 return(FALSE);
544 }
545 }
546 }
547 return(FALSE);
548 }
549 #define IS_ODD(j) ((j) & 1 )
550
551 double SIGN(double a,double b ) /* floating polong int SIGN transfer */
552 {
553 return( b < 0.0 ? -fabs(a) : fabs(a) );
554 }
555
556 double powi(double x, int n ) /* returns: x^n */
557 {
558 double p; /* holds partial product */
559
560 if( x == 0 )
561 return( 0. );
562
563 if( n < 0 ){ /* test for negative exponent */
564 n = -n;
565 x = 1/x;
566 }
567
568 p = IS_ODD(n) ? x : 1; /* test & set zero power */
569
570 while( n >>= 1 ){ /* now do the other powers */
571 x *= x; /* sq previous power of x */
572 if( IS_ODD(n) ) /* if low order bit set */
573 p *= x; /* then, multiply partial product by latest power of x */
574 }
575
576 return( p );
577 }
578