ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/kbond.c
Revision: 75
Committed: Mon Dec 15 18:03:29 2008 UTC (13 years, 5 months ago) by gilbertke
File size: 19542 byte(s)
Log Message:
fixes in typeing to pass MMFF validation suite
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 numeral(igeni, pa, 3);
145 numeral(igenk, pb, 3);
146 if( igeni <= igenk )
147 {
148 strcpy(pt_gen,pa);
149 strcat(pt_gen,pb);
150 }else
151 {
152 strcpy(pt_gen,pb);
153 strcat(pt_gen,pa);
154 }
155
156 ierr = FALSE;
157 // check 3 membered ring
158 if (rings.nring3 > 0 && bondk1.nbnd3 > 0 && ( is_ring31(ia) && is_ring31(ib) ) )
159 {
160 for (j=0; j < bondk1.nbnd3; j++)
161 {
162 if ( (strcmp(bondk1.kb3[j], pt) == 0) || (strcmp(bondk1.kb3[j],pt_gen) == 0) )
163 {
164 bonds_ff.bk[i] = bondk1.s3[j];
165 bonds_ff.bl[i] = bondk1.t3[j];
166 bonds_ff.index[i] = 0;
167 ierr = TRUE;
168 break;
169 }
170 }
171 if (ierr != TRUE)
172 {
173 Missing_constants = TRUE;
174 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],
175 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
176 fprintf(pcmlogfile,"%s\n",hatext);
177 //printf("%s\n",hatext);
178 }
179 }
180 if (ierr == TRUE)
181 goto L_10;
182 // check 4 membered ring
183 if (rings.nring4 > 0 && bondk1.nbnd4 > 0 && ( is_ring41(ia) && is_ring41(ib) ) )
184 {
185 for (j=0; j < bondk1.nbnd4; j++)
186 {
187 if ( (strcmp(bondk1.kb4[j],pt) == 0) || (strcmp(bondk1.kb4[j],pt_gen) == 0) )
188 {
189 bonds_ff.bk[i] = bondk1.s4[j];
190 bonds_ff.bl[i] = bondk1.t4[j];
191 bonds_ff.index[i] = 0;
192 ierr = TRUE;
193 break;
194 }
195 }
196 if (ierr != TRUE)
197 {
198 Missing_constants = TRUE;
199 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],
200 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
201 fprintf(pcmlogfile,"%s\n",hatext);
202 // printf("%s\n",hatext);
203 }
204 }
205 if (ierr == TRUE)
206 goto L_10;
207 // check 5 membered ring
208 if (rings.nring5 > 0 && bondk1.nbnd5 > 0 && (is_ring51(ia) != -1) && (is_ring51(ib) != -1))
209 {
210 for (j=0; j < bondk1.nbnd5; j++)
211 {
212 if (strcmp(bondk1.kb5[j],pt) == 0 )
213 {
214 bonds_ff.bk[i] = bondk1.s5[j];
215 bonds_ff.bl[i] = bondk1.t5[j];
216 bonds_ff.index[i] = 0;
217 ierr = TRUE;
218 break;
219 }
220 }
221 // don't fail on missing param - check for regular param first
222 }
223 if (ierr == TRUE)
224 goto L_10;
225 // delocalized bonds in MMFF94
226 if (field.type == MMFF94 && bondk1.ndeloc > 0)
227 {
228 if ( is_delocalbond(ia, ib) )
229 {
230 for (k = 0; k < MAXIAT; k++)
231 {
232 if (atom[ia].iat[k] == ib)
233 {
234 if (atom[ia].bo[k] == 1) // single bond
235 {
236 for (j=0; j < bondk1.ndeloc; j++)
237 {
238 if (strcmp(bondk1.kbdel[j],pt) == 0)
239 {
240 bonds_ff.bk[i] = bondk1.sdel[j];
241 bonds_ff.bl[i] = bondk1.tdel[j];
242 bonds_ff.index[i] = 1;
243 ierr = TRUE;
244 break;
245 }
246 }
247 }
248 }
249 }
250 }
251 if (ierr == TRUE)
252 goto L_10;
253 }
254 // regular parameters
255 for (j=0; j < bondk1.nbnd; j++)
256 {
257 if ( (strcmp(bondk1.kb[j], pt) == 0) )
258 {
259 bonds_ff.bk[i] = bondk1.s[j];
260 bonds_ff.bl[i] = bondk1.t[j];
261 bonds_ff.index[i] = 0;
262 ierr = TRUE;
263 if ( ((strcmp(pt," 2 2") == 0) || (strcmp(pt," 2 3") == 0))&& field.type == MMX)
264 {
265 for (k=0; k < MAXIAT; k++)
266 {
267 if (atom[ia].iat[k] == ib)
268 {
269 if (atom[ia].bo[k] == 1) // butadiene
270 {
271 bonds_ff.bl[i] = 1.48;
272 break;
273 }
274 }
275 }
276 }
277 break;
278 }
279 }
280 if (ierr == TRUE)
281 goto L_10;
282 //
283 // generalized parameters
284 for (j=0; j < bondk1.nbnd; j++)
285 {
286 if ( (strcmp(bondk1.kb[j], pt_gen) == 0) )
287 {
288 bonds_ff.bk[i] = bondk1.s[j];
289 bonds_ff.bl[i] = bondk1.t[j];
290 bonds_ff.index[i] = 0;
291 ierr = TRUE;
292 break;
293 }
294 }
295 if (ierr == TRUE)
296 goto L_10;
297 if (ierr != TRUE)
298 {
299 Missing_constants = TRUE;
300 sprintf(hatext,"Bond constants missing for bond %d - %d of type %d - %d\n",bonds_ff.i12[i][0],bonds_ff.i12[i][1],
301 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
302 fprintf(pcmlogfile,"%s\n",hatext);
303 // fprintf(errfile,"%s\n",hatext);
304 }
305 L_10:
306 continue;
307
308 }
309 }
310 if (field.type == MM3) // need check for electronegativity correction
311 {
312 for( i = 0; i < bonds_ff.nbnd; i++ )
313 {
314 ia = bonds_ff.i12[i][0];
315 ib = bonds_ff.i12[i][1];
316 iit = atom[bonds_ff.i12[i][0]].tclass;
317 kit = atom[bonds_ff.i12[i][1]].tclass;
318 if (iit > kit)
319 {
320 iend = ib;
321 kend = ia;
322 it = 1000*kit+iit;
323 }else
324 {
325 iend = ia;
326 kend = ib;
327 it = 1000*iit+kit;
328 }
329 numi = 0;
330 for (j=0; j < 6; j++)
331 {
332 iatt_electro[j] = 0;
333 iatt_const[j] = 0.0;
334 jatt_electro[j] = 0;
335 jatt_const[j] = 0.0;
336 }
337 for (j=0; j < electroneg.nelecti; j++)
338 {
339 if (it == electroneg.ibond[j])
340 {
341 for (k=0; k < MAXIAT; k++)
342 {
343 if (atom[iend].iat[k] != kend && atom[atom[iend].iat[k]].tclass == electroneg.iattach[j])
344 {
345 iatt_electro[numi]++;
346 iatt_const[numi] = electroneg.icorr[j];
347 numi++;
348 if (numi > 5) message_alert("error in kbond numi > 5","Error");
349 }
350 }
351 if (iit == kit)
352 {
353 for (k=0; k < MAXIAT; k++)
354 {
355 if (atom[kend].iat[k] != iend && atom[atom[kend].iat[k]].tclass == electroneg.iattach[j])
356 {
357 iatt_electro[numi]++;
358 iatt_const[numi] = electroneg.icorr[j];
359 numi++;
360 if (numi > 5) message_alert("error in kbond numi > 5","Error");
361 }
362 }
363 }
364 }
365 }
366 // loop over ib " " " "
367 numj = 0;
368 jatt_electro[0] = jatt_electro[1] = jatt_electro[2] = 0;
369 jatt_const[0] = jatt_const[1] = jatt_const[2] = 0.0;
370 for (j=0; j < electroneg.nelectj; j++)
371 {
372 if (it == electroneg.jbond[j])
373 {
374 for (k=0; k < MAXIAT; k++)
375 {
376 if (atom[kend].iat[k] != iend && atom[atom[kend].iat[k]].tclass == electroneg.jattach[j])
377 {
378 jatt_electro[numj]++;
379 jatt_const[numj] = electroneg.jcorr[j];
380 numj++;
381 if (numj > 5) message_alert("error in kbond numj > 5","Error");
382 }
383 }
384 }
385 }
386 // correct bond length
387 dbl = 0.0;
388 iy = -1;
389 for (j=0; j < numi; j++)
390 {
391 iy++;
392 dbl += iatt_const[j]*powi(0.62,iy);
393 }
394 for (j=0; j < numj; j++)
395 {
396 iy++;
397 dbl += jatt_const[j]*powi(0.62,iy);
398 }
399 bonds_ff.bl[i] += dbl;
400 //adjust force const
401 if (dbl != 0.0)
402 {
403 if (it == 1005 || it == 5022 || it == 5056)
404 {
405 bl0 = bonds_ff.bl[i] - dbl;
406 dks = deltaks(dbl,bl0);
407 bonds_ff.bk[i] += dks;
408 }
409 }
410
411 }
412 }
413 if (Missing_constants == TRUE)
414 {
415 fprintf(pcmlogfile,"Error assigning constants in: %s\n",Struct_Title);
416 fprintf(pcmlogfile,"%s\n",hatext);
417 return FALSE;
418 }
419
420 return TRUE;
421 }
422 // ==============================================
423 double deltaks(double dbl,double bl0)
424 {
425 double a,b,pi,clite,amu,conv,delta;
426 a = 1.3982;
427 b = -0.0001023;
428 pi = acos(-1.0);
429 clite = 2.9979e10;
430 amu = 1.660531e-24;
431 conv = 4.0*pi*pi*clite*clite*1.008*amu*1.0e-5;
432 delta = conv*dbl*2.0*(bl0-a)/(b*b);
433 return(delta);
434 }
435 /* ------------------------------------------ */
436 int find_bond(int ia, int ib)
437 {
438 int i;
439 int iz;
440
441 iz = -1;
442 for (i=0; i < bonds_ff.nbnd; i++)
443 {
444 if (ia == bonds_ff.i12[i][0] && ib == bonds_ff.i12[i][1])
445 return(i);
446 if (ib == bonds_ff.i12[i][0] && ia == bonds_ff.i12[i][1])
447 return(i);
448 }
449 return(iz);
450 }
451
452 /* ------------------------------------------ */
453 void angle(int k,int i,int imi,float *a13)
454 {
455 float disik, disim, dx1,dy1,dz1, dx2,dy2,dz2;
456
457 /* returns angle between k-i-imi */
458
459 dx1 = atom[i].x - atom[k].x;
460 dy1 = atom[i].y - atom[k].y;
461 dz1 = atom[i].z - atom[k].z;
462 dx2 = atom[i].x - atom[imi].x;
463 dy2 = atom[i].y - atom[imi].y;
464 dz2 = atom[i].z - atom[imi].z;
465
466 disik = sqrt( dx1*dx1 + dy1*dy1 + dz1*dz1);
467 disim = sqrt( dx2*dx2 + dy2*dy2 + dz2*dz2);
468
469 *a13 = angs( dx1, dy1, dz1, dx2, dy2, dz2, disik, disim );
470
471 *a13 *= 57.2958;
472
473 return;
474 }
475 /* ------------------------------ */
476 double angs(float x1,float y1,float z1,float x2,float y2,float z2,float d1,float d2)
477 {
478 double angs_v, cosa;
479
480 /* *** calculates angles for subroutine ebend *** */
481 cosa = (x1*x2 + y1*y2 + z1*z2)/(d1*d2);
482 if( fabs( cosa ) > 1.0 )
483 cosa /= fabs( cosa );
484 angs_v = arcos( cosa );
485 return( angs_v );
486 }
487 /* ----------------------- */
488 double arcos(float x)
489 {
490 double arcos_v, y;
491
492 y = x;
493 if( y > 1.0 )
494 y = 1.0;
495 if( y < -1.0 )
496 y = -1.0;
497 arcos_v = acos( y );
498 return( arcos_v );
499 }
500 /* ------------------------------------------ */
501 int is_delocalbond(int ia, int ib)
502 {
503 // test for delocalized bond
504 // if ia & ib are aromatic and part of 6 membered ring
505 // if not aromatic is bond order 1
506
507 long int aromatic_mask;
508 int i, j;
509 int jdbl, kdbl;
510
511 aromatic_mask = (1L << AROMATIC_MASK);
512
513 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
514 && (is_ring62(ia,ib) == TRUE) )
515 {
516 return(FALSE);
517 }
518
519 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
520 && (is_ring52(ia,ib) == TRUE) )
521 {
522 return(FALSE);
523 }
524 if (atom[ia].type == 37 && atom[ib].type == 37 && (is_ring62(ia,ib) != FALSE))
525 return FALSE;
526 if (atom[ia].type == 37 && atom[ib].type == 37 && (is_ring52(ia,ib) != FALSE))
527 return FALSE;
528
529 for (i=0; i < MAXIAT; i++)
530 {
531 if (atom[ia].iat[i] == ib)
532 {
533 if (atom[ia].bo[i] == 1)
534 {
535 jdbl = FALSE;
536 kdbl = FALSE;
537 for (j=0; j < MAXIAT; j++)
538 {
539 if (atom[ia].bo[j] >= 2 && atom[ia].bo[j] != 9)
540 jdbl = TRUE;
541 if (atom[ib].bo[j] >= 2 && atom[ib].bo[j] != 9)
542 kdbl = TRUE;
543 }
544 if (jdbl == TRUE && kdbl == TRUE)
545 return (TRUE);
546 else if (jdbl == TRUE && (atom[ib].flags & aromatic_mask) )
547 return (TRUE);
548 else if (kdbl == TRUE && (atom[ia].flags & aromatic_mask) )
549 return (TRUE);
550 else
551 return(FALSE);
552 }
553 }
554 }
555 return(FALSE);
556 }
557 #define IS_ODD(j) ((j) & 1 )
558
559 double SIGN(double a,double b ) /* floating polong int SIGN transfer */
560 {
561 return( b < 0.0 ? -fabs(a) : fabs(a) );
562 }
563
564 double powi(double x, int n ) /* returns: x^n */
565 {
566 double p; /* holds partial product */
567
568 if( x == 0 )
569 return( 0. );
570
571 if( n < 0 ){ /* test for negative exponent */
572 n = -n;
573 x = 1/x;
574 }
575
576 p = IS_ODD(n) ? x : 1; /* test & set zero power */
577
578 while( n >>= 1 ){ /* now do the other powers */
579 x *= x; /* sq previous power of x */
580 if( IS_ODD(n) ) /* if low order bit set */
581 p *= x; /* then, multiply partial product by latest power of x */
582 }
583
584 return( p );
585 }
586