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 (12 years, 4 months ago) by tjod
File size: 24046 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #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