ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/kbond.c
Revision: 65
Committed: Fri Dec 5 02:21:07 2008 UTC (12 years, 9 months ago) by gilbertke
File size: 19168 byte(s)
Log Message:
added call to type in rd_sdf, fixed problem in type_mmx
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 L_10:
298 continue;
299
300 }
301 }
302 if (field.type == MM3) // need check for electronegativity correction
303 {
304 for( i = 0; i < bonds_ff.nbnd; i++ )
305 {
306 ia = bonds_ff.i12[i][0];
307 ib = bonds_ff.i12[i][1];
308 iit = atom[bonds_ff.i12[i][0]].tclass;
309 kit = atom[bonds_ff.i12[i][1]].tclass;
310 if (iit > kit)
311 {
312 iend = ib;
313 kend = ia;
314 it = 1000*kit+iit;
315 }else
316 {
317 iend = ia;
318 kend = ib;
319 it = 1000*iit+kit;
320 }
321 numi = 0;
322 for (j=0; j < 6; j++)
323 {
324 iatt_electro[j] = 0;
325 iatt_const[j] = 0.0;
326 jatt_electro[j] = 0;
327 jatt_const[j] = 0.0;
328 }
329 for (j=0; j < electroneg.nelecti; j++)
330 {
331 if (it == electroneg.ibond[j])
332 {
333 for (k=0; k < MAXIAT; k++)
334 {
335 if (atom[iend].iat[k] != kend && atom[atom[iend].iat[k]].tclass == electroneg.iattach[j])
336 {
337 iatt_electro[numi]++;
338 iatt_const[numi] = electroneg.icorr[j];
339 numi++;
340 if (numi > 5) message_alert("error in kbond numi > 5","Error");
341 }
342 }
343 if (iit == kit)
344 {
345 for (k=0; k < MAXIAT; k++)
346 {
347 if (atom[kend].iat[k] != iend && atom[atom[kend].iat[k]].tclass == electroneg.iattach[j])
348 {
349 iatt_electro[numi]++;
350 iatt_const[numi] = electroneg.icorr[j];
351 numi++;
352 if (numi > 5) message_alert("error in kbond numi > 5","Error");
353 }
354 }
355 }
356 }
357 }
358 // loop over ib " " " "
359 numj = 0;
360 jatt_electro[0] = jatt_electro[1] = jatt_electro[2] = 0;
361 jatt_const[0] = jatt_const[1] = jatt_const[2] = 0.0;
362 for (j=0; j < electroneg.nelectj; j++)
363 {
364 if (it == electroneg.jbond[j])
365 {
366 for (k=0; k < MAXIAT; k++)
367 {
368 if (atom[kend].iat[k] != iend && atom[atom[kend].iat[k]].tclass == electroneg.jattach[j])
369 {
370 jatt_electro[numj]++;
371 jatt_const[numj] = electroneg.jcorr[j];
372 numj++;
373 if (numj > 5) message_alert("error in kbond numj > 5","Error");
374 }
375 }
376 }
377 }
378 // correct bond length
379 dbl = 0.0;
380 iy = -1;
381 for (j=0; j < numi; j++)
382 {
383 iy++;
384 dbl += iatt_const[j]*powi(0.62,iy);
385 }
386 for (j=0; j < numj; j++)
387 {
388 iy++;
389 dbl += jatt_const[j]*powi(0.62,iy);
390 }
391 bonds_ff.bl[i] += dbl;
392 //adjust force const
393 if (dbl != 0.0)
394 {
395 if (it == 1005 || it == 5022 || it == 5056)
396 {
397 bl0 = bonds_ff.bl[i] - dbl;
398 dks = deltaks(dbl,bl0);
399 bonds_ff.bk[i] += dks;
400 }
401 }
402
403 }
404 }
405 if (Missing_constants == TRUE)
406 {
407 fprintf(pcmlogfile,"Error assigning constants in: %s\n",Struct_Title);
408 fprintf(pcmlogfile,"%s\n",hatext);
409 return FALSE;
410 }
411
412 return TRUE;
413 }
414 // ==============================================
415 double deltaks(double dbl,double bl0)
416 {
417 double a,b,pi,clite,amu,conv,delta;
418 a = 1.3982;
419 b = -0.0001023;
420 pi = acos(-1.0);
421 clite = 2.9979e10;
422 amu = 1.660531e-24;
423 conv = 4.0*pi*pi*clite*clite*1.008*amu*1.0e-5;
424 delta = conv*dbl*2.0*(bl0-a)/(b*b);
425 return(delta);
426 }
427 /* ------------------------------------------ */
428 int find_bond(int ia, int ib)
429 {
430 int i;
431 int iz;
432
433 iz = -1;
434 for (i=0; i < bonds_ff.nbnd; i++)
435 {
436 if (ia == bonds_ff.i12[i][0] && ib == bonds_ff.i12[i][1])
437 return(i);
438 if (ib == bonds_ff.i12[i][0] && ia == bonds_ff.i12[i][1])
439 return(i);
440 }
441 return(iz);
442 }
443
444 /* ------------------------------------------ */
445 void angle(int k,int i,int imi,float *a13)
446 {
447 float disik, disim, dx1,dy1,dz1, dx2,dy2,dz2;
448
449 /* returns angle between k-i-imi */
450
451 dx1 = atom[i].x - atom[k].x;
452 dy1 = atom[i].y - atom[k].y;
453 dz1 = atom[i].z - atom[k].z;
454 dx2 = atom[i].x - atom[imi].x;
455 dy2 = atom[i].y - atom[imi].y;
456 dz2 = atom[i].z - atom[imi].z;
457
458 disik = sqrt( dx1*dx1 + dy1*dy1 + dz1*dz1);
459 disim = sqrt( dx2*dx2 + dy2*dy2 + dz2*dz2);
460
461 *a13 = angs( dx1, dy1, dz1, dx2, dy2, dz2, disik, disim );
462
463 *a13 *= 57.2958;
464
465 return;
466 }
467 /* ------------------------------ */
468 double angs(float x1,float y1,float z1,float x2,float y2,float z2,float d1,float d2)
469 {
470 double angs_v, cosa;
471
472 /* *** calculates angles for subroutine ebend *** */
473 cosa = (x1*x2 + y1*y2 + z1*z2)/(d1*d2);
474 if( fabs( cosa ) > 1.0 )
475 cosa /= fabs( cosa );
476 angs_v = arcos( cosa );
477 return( angs_v );
478 }
479 /* ----------------------- */
480 double arcos(float x)
481 {
482 double arcos_v, y;
483
484 y = x;
485 if( y > 1.0 )
486 y = 1.0;
487 if( y < -1.0 )
488 y = -1.0;
489 arcos_v = acos( y );
490 return( arcos_v );
491 }
492 /* ------------------------------------------ */
493 int is_delocalbond(int ia, int ib)
494 {
495 // test for delocalized bond
496 // if ia & ib are aromatic and part of 6 membered ring
497 // if not aromatic is bond order 1
498
499 long int aromatic_mask;
500 int i, j;
501 int jdbl, kdbl;
502
503 aromatic_mask = (1L << AROMATIC_MASK);
504
505 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
506 && (is_ring62(ia,ib) == TRUE) )
507 {
508 return(FALSE);
509 }
510
511 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
512 && (is_ring52(ia,ib) == TRUE) )
513 {
514 return(FALSE);
515 }
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