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