ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/kbond.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 24046 byte(s)
Log Message:
test

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 #include "fix.h"
9
10 FILE * fopen_path ( char * , char * , char * ) ;
11 void find_metalbond(int *, int, int,char *,float *, float *);
12 void angle(int ,int ,int ,float *);
13 int is_ring31(int);
14 int is_ring41(int);
15 int is_ring51(int);
16 int is_ring61(int);
17 int is_ring42(int, int);
18 int is_ring52(int, int);
19 int is_ring62(int,int);
20 void pcmfout(int);
21 double arcos(float );
22 double angs(float,float,float,float,float ,float ,float ,float);
23 double powi(double, int);
24 void numeral(int, char *, int);
25 int is_delocalbond(int ia, int ib);
26 void message_alert(char *, char *);
27 double sign(double ,double );
28 void transbond(int *,char *,float *,float *);
29 int is_bond(int,int);
30 double deltaks(double,double);
31 float get_bond(int, int);
32
33
34
35 EXTERN struct t_bondk1 {
36 int use_ring3, use_ring4, use_ring5;
37 int nbnd, nbnd3, nbnd4, nbnd5, ndeloc;
38 char kb[MAXBONDCONST][7], kb3[MAXBOND3CONST][7],kb4[MAXBOND4CONST][7],kb5[MAXBOND5CONST][7];
39 char kbdel[MAXBONDDELOC][7];
40 float s[MAXBONDCONST], t[MAXBONDCONST];
41 float s3[MAXBOND3CONST], t3[MAXBOND3CONST];
42 float s4[MAXBOND4CONST], t4[MAXBOND4CONST];
43 float s5[MAXBOND5CONST], t5[MAXBOND5CONST];
44 float sdel[MAXBONDDELOC], tdel[MAXBONDDELOC];
45 } bondk1;
46
47
48 EXTERN struct t_bondpk {
49 int npibond,npibond44i,npibond45i,npibond46i,npibond55i,npibond56i,npibond66i;
50 int npibond44o,npibond45o,npibond46o,npibond55o,npibond56o,npibond66o;
51 int npibond40,npibond50,npibond60;
52 int use_pibond,use_pibond44i,use_pibond45i,use_pibond46i,use_pibond55i,use_pibond56i,use_pibond66i;
53 int use_pibond44o,use_pibond45o,use_pibond46o,use_pibond55o,use_pibond56o,use_pibond66o;
54 int use_pibond40,use_pibond50,use_pibond60;
55 char kb[MAXPIBONDCONST][7],kb44i[20][7],kb45i[20][7],kb46i[20][7],kb55i[20][7],kb56i[20][7],kb66i[20][7];
56 char kb44o[20][7],kb45o[20][7],kb46o[20][7],kb55o[20][7],kb56o[20][7],kb66o[20][7];
57 char kb40[20][7],kb50[20][7],kb60[20][7];
58 float bk[MAXPIBONDCONST], bl[MAXPIBONDCONST], bmom[MAXPIBONDCONST],sslop[MAXPIBONDCONST], tslop[MAXPIBONDCONST], tslop2[MAXPIBONDCONST];
59 float bk44i[20], bl44i[20], bmom44i[20],sslop44i[20], tslop44i[20], tslop244i[20];
60 float bk45i[20], bl45i[20], bmom45i[20],sslop45i[20], tslop45i[20], tslop245i[20];
61 float bk46i[20], bl46i[20], bmom46i[20],sslop46i[20], tslop46i[20], tslop246i[20];
62 float bk55i[20], bl55i[20], bmom55i[20],sslop55i[20], tslop55i[20], tslop255i[20];
63 float bk56i[20], bl56i[20], bmom56i[20],sslop56i[20], tslop56i[20], tslop256i[20];
64 float bk66i[20], bl66i[20], bmom66i[20],sslop66i[20], tslop66i[20], tslop266i[20];
65 float bk44o[20], bl44o[20], bmom44o[20],sslop44o[20], tslop44o[20], tslop244o[20];
66 float bk45o[20], bl45o[20], bmom45o[20],sslop45o[20], tslop45o[20], tslop245o[20];
67 float bk46o[20], bl46o[20], bmom46o[20],sslop46o[20], tslop46o[20], tslop246o[20];
68 float bk55o[20], bl55o[20], bmom55o[20],sslop55o[20], tslop55o[20], tslop255o[20];
69 float bk56o[20], bl56o[20], bmom56o[20],sslop56o[20], tslop56o[20], tslop256o[20];
70 float bk66o[20], bl66o[20], bmom66o[20],sslop66o[20], tslop66o[20], tslop266o[20];
71 float bk40[20], bl40[20], bmom40[20],sslop40[20], tslop40[20], tslop240[20];
72 float bk50[20], bl50[20], bmom50[20],sslop50[20], tslop50[20], tslop250[20];
73 float bk60[20], bl60[20], bmom60[20],sslop60[20], tslop60[20], tslop260[20];
74 } bondpk;
75
76
77 EXTERN struct t_electroneg {
78 int nelecti, nelectj;
79 int itype[200],ibond[200],iattach[200];
80 int jtype[200],jbond[200],jattach[200];
81 float icorr[200],jcorr[200];
82 } electroneg;
83
84 EXTERN struct t_ts_bondorder {
85 float fbnd[15];
86 } ts_bondorder;
87 EXTERN struct t_files {
88 int nfiles, append, batch, icurrent;
89 int istrselect;
90 } files;
91
92 EXTERN int Missing_constants;
93 //EXTERN FILE *errfile;
94
95 int kbond()
96 {
97 long int pi_mask;
98 int ia4,ia5,ia6, ib4,ib5,ib6, iab4, iab5, iab6;
99 int i, j,k, iit, kit, it, ierr,nbpi, ii, kk, i1, i2;
100 int igeni, igenk, ifound;
101 int ia, ib, iend, kend, iy;
102 int numi, numj, iatt_electro[6], jatt_electro[6];
103 float jatt_const[6],iatt_const[6];
104 float bconst, blength;
105 float sslope, tslope, tslope2;
106 double dbl,dks,bl0;
107 char hatext[80];
108 char pa[4],pb[4], pt[7], pt_gen[7];
109
110 pi_mask = 1L << PI_MASK; /* this mask is for pi atoms - now uses bit 0 */
111 nbpi = 0;
112 if( natom != 0 )
113 {
114
115 for( i = 0; i < bonds_ff.nbnd; i++ )
116 {
117 ia = bonds_ff.i12[i][0];
118 ib = bonds_ff.i12[i][1];
119 iit = atom[bonds_ff.i12[i][0]].type;
120 kit = atom[bonds_ff.i12[i][1]].type;
121
122 igeni = atom[bonds_ff.i12[i][0]].tclass;
123 igenk = atom[bonds_ff.i12[i][1]].tclass;
124
125 // metal general class
126 if (iit > 300)
127 igeni = 300;
128 if (kit > 300)
129 igenk = 300;
130
131 if( field.type == MMX)
132 {
133 if (iit == 40 && ( kit != 2 && kit != 3 && kit != 4 && kit != 40) )
134 iit = 2;
135 if (kit == 40 && ( iit != 2 && iit != 3 && iit != 4 && iit != 40) )
136 kit = 2;
137 }
138
139 numeral(iit, pa, 3);
140 numeral(kit, pb, 3);
141 strcpy(pt,pa);
142 strcat(pt,pb);
143 if( iit <= kit )
144 {
145 strcpy(pt,pa);
146 strcat(pt,pb);
147 }else
148 {
149 strcpy(pt,pb);
150 strcat(pt,pa);
151 }
152
153 numeral(igeni, pa, 3);
154 numeral(igenk, pb, 3);
155 if( igeni <= igenk )
156 {
157 strcpy(pt_gen,pa);
158 strcat(pt_gen,pb);
159 }else
160 {
161 strcpy(pt_gen,pb);
162 strcat(pt_gen,pa);
163 }
164
165 ierr = FALSE;
166 // check 3 membered ring
167 if (rings.nring3 > 0 && bondk1.nbnd3 > 0 && ( is_ring31(ia) && is_ring31(ib) ) )
168 {
169 for (j=0; j < bondk1.nbnd3; j++)
170 {
171 if ( (strcmp(bondk1.kb3[j], pt) == 0) || (strcmp(bondk1.kb3[j],pt_gen) == 0) )
172 {
173 bonds_ff.bk[i] = bondk1.s3[j];
174 bonds_ff.bl[i] = bondk1.t3[j];
175 bonds_ff.index[i] = 0;
176 ierr = TRUE;
177 break;
178 }
179 }
180 if (ierr != TRUE)
181 {
182 Missing_constants = TRUE;
183 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],
184 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
185 }
186 }
187 if (ierr == TRUE)
188 goto L_10;
189 // check 4 membered ring
190 if (rings.nring4 > 0 && bondk1.nbnd4 > 0 && ( is_ring41(ia) && is_ring41(ib) ) )
191 {
192 for (j=0; j < bondk1.nbnd4; j++)
193 {
194 if ( (strcmp(bondk1.kb4[j],pt) == 0) || (strcmp(bondk1.kb4[j],pt_gen) == 0) )
195 {
196 bonds_ff.bk[i] = bondk1.s4[j];
197 bonds_ff.bl[i] = bondk1.t4[j];
198 bonds_ff.index[i] = 0;
199 ierr = TRUE;
200 break;
201 }
202 }
203 if (ierr != TRUE)
204 {
205 Missing_constants = TRUE;
206 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],
207 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
208 }
209 }
210 if (ierr == TRUE)
211 goto L_10;
212 // check 5 membered ring
213 if (rings.nring5 > 0 && bondk1.nbnd5 > 0 && (is_ring51(ia) != -1) && (is_ring51(ib) != -1))
214 {
215 for (j=0; j < bondk1.nbnd5; j++)
216 {
217 if (strcmp(bondk1.kb5[j],pt) == 0 )
218 {
219 bonds_ff.bk[i] = bondk1.s5[j];
220 bonds_ff.bl[i] = bondk1.t5[j];
221 bonds_ff.index[i] = 0;
222 ierr = TRUE;
223 break;
224 }
225 }
226 // don't fail on missing param - check for regular param first
227 }
228 if (ierr == TRUE)
229 goto L_10;
230 // delocalized bonds in MMFF94
231 if (field.type == MMFF94 && bondk1.ndeloc > 0)
232 {
233 if ( is_delocalbond(ia, ib) )
234 {
235 for (k = 0; k < MAXIAT; k++)
236 {
237 if (atom[ia].iat[k] == ib)
238 {
239 if (atom[ia].bo[k] == 1) // single bond
240 {
241 for (j=0; j < bondk1.ndeloc; j++)
242 {
243 if (strcmp(bondk1.kbdel[j],pt) == 0)
244 {
245 bonds_ff.bk[i] = bondk1.sdel[j];
246 bonds_ff.bl[i] = bondk1.tdel[j];
247 bonds_ff.index[i] = 1;
248 ierr = TRUE;
249 break;
250 }
251 }
252 }
253 }
254 }
255 }
256 if (ierr == TRUE)
257 goto L_10;
258 }
259 // regular parameters
260 for (j=0; j < bondk1.nbnd; j++)
261 {
262 if ( (strcmp(bondk1.kb[j], pt) == 0) )
263 {
264 bonds_ff.bk[i] = bondk1.s[j];
265 bonds_ff.bl[i] = bondk1.t[j];
266 bonds_ff.index[i] = 0;
267 ierr = TRUE;
268 if ( ((strcmp(pt," 2 2") == 0) || (strcmp(pt," 2 3") == 0))&& field.type == MMX)
269 {
270 for (k=0; k < MAXIAT; k++)
271 {
272 if (atom[ia].iat[k] == ib)
273 {
274 if (atom[ia].bo[k] == 1) // butadiene
275 {
276 bonds_ff.bl[i] = 1.48;
277 break;
278 }
279 }
280 }
281 }
282 break;
283 }
284 }
285 if (ierr == TRUE)
286 goto L_10;
287 //
288 // generalized parameters
289 for (j=0; j < bondk1.nbnd; j++)
290 {
291 if ( (strcmp(bondk1.kb[j], pt_gen) == 0) )
292 {
293 bonds_ff.bk[i] = bondk1.s[j];
294 bonds_ff.bl[i] = bondk1.t[j];
295 bonds_ff.index[i] = 0;
296 ierr = TRUE;
297 break;
298 }
299 }
300 if (ierr == TRUE)
301 goto L_10;
302 // check for metal bonds - have not found specific metal bonds - now looking for general mmx type bonds
303 if ( iit >= 300 || kit >= 300 )
304 {
305 find_metalbond(&ierr,bonds_ff.i12[i][0],bonds_ff.i12[i][1], pt_gen, &bconst, &blength);
306 if (ierr == TRUE)
307 {
308 bonds_ff.bk[i] = bconst;
309 bonds_ff.bl[i] = blength;
310 }
311 }
312 if (ierr == TRUE)
313 {
314 fprintf(pcmoutfile,"Using generalized metal bonds for MMX force field - %d %d of types %d %d\n",ia,ib,iit,kit);
315 goto L_10;
316 }
317 if (ierr != TRUE)
318 {
319 Missing_constants = TRUE;
320 sprintf(hatext,"Bond constants missing for bond %d - %d of type %d - %d\n",bonds_ff.i12[i][0],bonds_ff.i12[i][1],
321 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
322 }
323 L_10:
324 continue;
325
326 }
327 }
328 if (Missing_constants == TRUE)
329 {
330 return FALSE;
331 }
332
333 if (fixdis.FixBond) // fixed bonds replace bond constant with fixed constant
334 {
335 for (i=0; i < fixdis.nfxstr; i++)
336 {
337 ia = fixdis.ifxstr[i][0];
338 ib = fixdis.ifxstr[i][1];
339 if (is_bond(ia,ib))
340 {
341 for (j=0; j < bonds_ff.nbnd; j++)
342 {
343 i1 = bonds_ff.i12[j][0];
344 i2 = bonds_ff.i12[j][1];
345 if ( ((ia == i1) && (ib == i2)) || ((ia == i2) && (ib == i1)) )
346 {
347 bonds_ff.bl[j] = fixdis.fxstrd[i];
348 bonds_ff.bk[j] = fixdis.fxstrc[i];
349 }
350 }
351 }
352 }
353 }
354
355 return TRUE;
356 }
357 /* ------------------------------------------ */
358 void find_metalbond(int *ierr,int at1, int at2, char *pt, float *bconst, float *blength)
359 {
360 int j, k, itm, ito, nusi, ntm, mi,imi;
361 long int mask5, mask6, mask7, mask8;
362 float a13;
363
364 *blength = 0.0;
365
366 for (j=0; j < bondk1.nbnd; j++)
367 {
368 if (strcmp(bondk1.kb[j], pt) == 0)
369 {
370 *bconst = bondk1.s[j];
371 *blength = bondk1.t[j];
372 *ierr = TRUE;
373 break;
374 }
375 }
376 if (strcmp(pt,"300300") == 0)
377 {
378 if (atom[at1].type == atom[at2].type)
379 {
380 *blength += 2*atom[at1].radius - 2.52;
381 } else
382 {
383 *blength += atom[at1].radius + atom[at2].radius -2.52;
384 }
385 for (k=0; k < MAXIAT; k++)
386 {
387 if (atom[at1].iat[k] == at2)
388 {
389 if (atom[at1].bo[k] > 1 && atom[at1].bo[k] != 9)
390 *blength += -atom[at1].bo[k]*0.19;
391 }
392 }
393 }
394 if ( (atom[at1].type >= 300 && atom[at2].type < 300 ) || (atom[at2].type >= 300 && atom[at1].type < 300) )
395 {
396 if (atom[at1].type >= 300)
397 {
398 itm = at1;
399 ito = at2;
400 }else
401 {
402 itm = at2;
403 ito = at1;
404 }
405 if (*blength == 0.00)
406 *blength += atom[itm].vdw_radius - 1.26;
407 /* flags 5 and 6 are used to mark electron count
408 * 1,0 for saturated: 0,1 gt 18 electron and 1,1 for unsaturated
409 *
410 * flags 7 and 8 are used for high spin, low spin and sq planar
411 * 1,0 for low spin: 0,1 for sq planar : 1,1, for high spin
412 *
413 * find out if unsaturated and then the spin state
414 *
415 * nusi - 0 : coord sat'd (also high spin)
416 * nusi - 1 : >18 electrons
417 * nusi - 2 : low spin
418 * nusi - 3 : square planar */
419 mask5 = 1L << SATMET_MASK;
420 mask6 = 1L << GT18e_MASK;
421 mask7 = 1L << LOWSPIN_MASK;
422 mask8 = 1L << SQPLAN_MASK;
423 nusi = 0;
424 if( ( atom[itm].flags & mask6 ) && !( atom[itm].flags & mask5 ) )
425 {
426 nusi = 1;
427 }else if( ( atom[itm].flags & mask5 ) && ( atom[itm].flags & mask6 ) )
428 {
429 if( ( atom[itm].flags & mask7 ) && !( atom[itm].flags & mask8 ) )
430 {
431 nusi = 2;
432 }else if(( atom[itm].flags & mask8 ) && !( atom[itm].flags & mask7 ) )
433 {
434 nusi = 3;
435 }
436 }
437
438 if( nusi == 2 || nusi == 3 )
439 {
440 /* unsaturated carbon */
441 if( atom[ito].type == 2 )
442 {
443 *blength -= 0.15;
444 /* unsaturated oxygen or sulfur */
445 }else if( atom[ito].type == 6 || atom[ito].type == 15 )
446 {
447 *blength -= 0.1;
448 /* unsaturated halogens */
449 }else if( atom[ito].type >= 11 && atom[ito].type <= 15 )
450 {
451 *blength -= 0.2;
452 }
453 *bconst *= 1.5;
454 }
455 if( nusi == 1 )
456 *blength += .15;
457
458 /* sq planar complexes
459 * found a metal-halogen search other bonds to metal for trans effects */
460 if( nusi == 3 && (atom[ito].type== 12 || atom[ito].type == 13) )
461 {
462 for( mi = 0; mi < MAXIAT; mi++ )
463 {
464 imi = atom[itm].iat[mi];
465 if( atom[itm].iat[mi] != 0 )
466 {
467 ntm = atom[imi].type;
468 if( ntm == 1 || ntm == 2 || ntm == 4 || ntm == 5 || ntm == 25 || ntm == 30 )
469 {
470 angle( ito, itm, imi,&a13 );
471 if( a13 > 150 || a13 < -150. )
472 {
473 *blength += .20;
474 fprintf(pcmoutfile,"Trans influence (+0.20 a) on metal-cl or br bond for %d - %d\n",itm,ito);
475 }
476 }
477 }
478 }
479 }
480
481 /* more sq planar tests */
482 if( nusi == 3 && (atom[ito].type == 1 || atom[ito].type == 2 || atom[ito].type == 5 ||
483 atom[ito].type == 25 || atom[ito].type == 30) )
484 {
485 for( mi = 0; mi < MAXIAT; mi++ )
486 {
487 imi = atom[itm].iat[mi];
488 if( imi != 0 )
489 {
490 ntm = atom[imi].type;
491 if( ntm == 12 || ntm == 13 )
492 {
493 angle( ito, itm, imi,&a13);
494 if( a13 > 150 || a13 < -150. )
495 {
496 *blength -= .1;
497 fprintf(pcmoutfile,"trans influence (-0.1 a) on metal-ligand bond: %d - %d\n",itm,ito);
498 }
499 }
500 }
501 }
502 }
503 /* metal to nitrile */
504 if( atom[ito].type == 10 )
505 {
506 for( mi = 0; mi < MAXIAT; mi++ )
507 {
508 if( atom[ito].iat[mi] != 0 )
509 {
510 ntm = atom[atom[ito].iat[mi]].type;
511 if( ntm == 4 || ntm == 10 )
512 {
513 *blength = 1.91;
514 *bconst = 2.0;
515 }
516 }
517 }
518 }
519 /* metal to amine */
520 if( atom[ito].type == 7 )
521 {
522 for( mi = 0; mi < MAXIAT; mi++ )
523 {
524 if( atom[ito].iat[mi] != 0 )
525 {
526 if( atom[atom[ito].iat[mi]].type == 3 )
527 {
528 *blength = 1.70;
529 *bconst = 2.0;
530 }
531 }
532 }
533 }
534 /* metal to ammonium */
535 if( atom[ito].type == 41 )
536 {
537 for( mi = 0; mi < MAXIAT; mi++ )
538 {
539 if( atom[ito].iat[mi] != 0 )
540 {
541 ntm = atom[atom[ito].iat[mi]].type;
542 if( ntm == 7 )
543 *blength = 1.67;
544 if( ntm == 37 || ntm == 41 )
545 *blength = 1.73;
546 }
547 }
548 }
549 }
550
551 }
552 /* ------------------------------------------ */
553 int find_bond(int ia, int ib)
554 {
555 int i;
556 int iz;
557
558 iz = -1;
559 for (i=0; i < bonds_ff.nbnd; i++)
560 {
561 if (ia == bonds_ff.i12[i][0] && ib == bonds_ff.i12[i][1])
562 return(i);
563 if (ib == bonds_ff.i12[i][0] && ia == bonds_ff.i12[i][1])
564 return(i);
565 }
566 return(iz);
567 }
568
569 /* ------------------------------------------ */
570 void angle(int k,int i,int imi,float *a13)
571 {
572 float disik, disim, dx1,dy1,dz1, dx2,dy2,dz2;
573
574 /* returns angle between k-i-imi */
575
576 dx1 = atom[i].x - atom[k].x;
577 dy1 = atom[i].y - atom[k].y;
578 dz1 = atom[i].z - atom[k].z;
579 dx2 = atom[i].x - atom[imi].x;
580 dy2 = atom[i].y - atom[imi].y;
581 dz2 = atom[i].z - atom[imi].z;
582
583 disik = sqrt( dx1*dx1 + dy1*dy1 + dz1*dz1);
584 disim = sqrt( dx2*dx2 + dy2*dy2 + dz2*dz2);
585
586 *a13 = angs( dx1, dy1, dz1, dx2, dy2, dz2, disik, disim );
587
588 *a13 *= 57.2958;
589
590 return;
591 }
592 /* ------------------------------ */
593 double angs(float x1,float y1,float z1,float x2,float y2,float z2,float d1,float d2)
594 {
595 double angs_v, cosa;
596
597 /* *** calculates angles for subroutine ebend *** */
598 cosa = (x1*x2 + y1*y2 + z1*z2)/(d1*d2);
599 if( fabs( cosa ) > 1.0 )
600 cosa /= fabs( cosa );
601 angs_v = arcos( cosa );
602 return( angs_v );
603 }
604 /* ----------------------- */
605 double arcos(float x)
606 {
607 double arcos_v, y;
608
609 y = x;
610 if( y > 1.0 )
611 y = 1.0;
612 if( y < -1.0 )
613 y = -1.0;
614 arcos_v = acos( y );
615 return( arcos_v );
616 }
617 /* ------------------------------------------ */
618 int is_delocalbond(int ia, int ib)
619 {
620 // test for delocalized bond
621 // if ia & ib are aromatic and part of 6 membered ring
622 // if not aromatic is bond order 1
623
624 long int aromatic_mask;
625 int i, j;
626 int jdbl, kdbl;
627
628 aromatic_mask = (1L << AROMATIC_MASK);
629
630 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
631 && (is_ring62(ia,ib) != FALSE) )
632 return(FALSE);
633
634 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
635 && (is_ring52(ia,ib) != FALSE) )
636 return(FALSE);
637
638 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask) )
639 return(TRUE);
640 if (atom[ia].type == 37 && atom[ib].type == 37 && (is_ring62(ia,ib) != FALSE))
641 return FALSE;
642 if (atom[ia].type == 37 && atom[ib].type == 37 && (is_ring52(ia,ib) != FALSE))
643 return FALSE;
644
645 for (i=0; i < MAXIAT; i++)
646 {
647 if (atom[ia].iat[i] == ib)
648 {
649 if (atom[ia].bo[i] == 1)
650 {
651 jdbl = FALSE;
652 kdbl = FALSE;
653 for (j=0; j < MAXIAT; j++)
654 {
655 if (atom[ia].bo[j] >= 2 && atom[ia].bo[j] != 9)
656 jdbl = TRUE;
657 if (atom[ib].bo[j] >= 2 && atom[ib].bo[j] != 9)
658 kdbl = TRUE;
659 }
660 if (jdbl == TRUE && kdbl == TRUE)
661 return (TRUE);
662 else if (jdbl == TRUE && (atom[ib].flags & aromatic_mask) )
663 return (TRUE);
664 else if (kdbl == TRUE && (atom[ia].flags & aromatic_mask) )
665 return (TRUE);
666 else
667 return(FALSE);
668 }
669 }
670 }
671 return(FALSE);
672 }
673 #define IS_ODD(j) ((j) & 1 )
674
675 double sign(double a,double b ) /* floating polong int sign transfer */
676 {
677 return( b < 0.0 ? -fabs(a) : fabs(a) );
678 }
679
680 double powi(double x, int n ) /* returns: x^n */
681 {
682 double p; /* holds partial product */
683
684 if( x == 0 )
685 return( 0. );
686
687 if( n < 0 ){ /* test for negative exponent */
688 n = -n;
689 x = 1/x;
690 }
691
692 p = IS_ODD(n) ? x : 1; /* test & set zero power */
693
694 while( n >>= 1 ){ /* now do the other powers */
695 x *= x; /* sq previous power of x */
696 if( IS_ODD(n) ) /* if low order bit set */
697 p *= x; /* then, multiply partial product by latest power of x */
698 }
699
700 return( p );
701 }
702 /* ================================ */
703 float get_bond(int ia, int ib)
704 {
705 int i, ita,itb;
706 char pa[4],pb[4], pt[7];
707
708 ita = atom[ia].type;
709 itb = atom[ib].type;
710 numeral(ita,pa,3);
711 numeral(itb,pb,3);
712 if (ita <= itb)
713 {
714 strcpy(pt,pa);
715 strcat(pt,pb);
716 } else
717 {
718 strcpy(pt,pb);
719 strcat(pt,pa);
720 }
721 for (i=0; i < bondk1.nbnd; i++)
722 {
723 if ( (strcmp(bondk1.kb[i],pt) == 0))
724 {
725 return (bondk1.t[i]);
726 }
727 }
728 return 1.53; // default to c-c single
729 }
730