ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/kbond.c
(Generate patch)
# Line 6 | Line 6
6   #include "rings.h"
7   #include "field.h"
8  
9 FILE * fopen_path ( char * , char * , char * ) ;
10 void angle(int ,int ,int ,float *);
9   int is_ring31(int);
10   int is_ring41(int);
11   int is_ring51(int);
# Line 15 | Line 13
13   int is_ring42(int, int);
14   int is_ring52(int, int);
15   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);
16   void numeral(int, char *, int);
17   int is_delocalbond(int ia, int ib);
18   void message_alert(char *, char *);
19   double SIGN(double ,double  );
20   int is_bond(int,int);
27 double deltaks(double,double);
21   float get_bond(int, int);
22  
23  
# Line 39 | Line 32
32          float s5[MAXBOND5CONST], t5[MAXBOND5CONST];
33          float sdel[MAXBONDDELOC], tdel[MAXBONDDELOC];
34          }  bondk1;
35 <                
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 <        
35 >                        
36   EXTERN int Missing_constants;
85 EXTERN FILE *errfile;
37  
38   int kbond()
39   {
# Line 91 | Line 42
42      int i, j,k, iit, kit, it, ierr,nbpi, ii, kk, i1, i2;
43      int igeni, igenk, ifound;
44      int ia, ib, iend, kend, iy;
45 <    int numi, numj, iatt_electro[6], jatt_electro[6];
95 <    float jatt_const[6],iatt_const[6];
45 >    int numi, numj;
46      float bconst, blength;
47      float sslope, tslope, tslope2;
48      double dbl,dks,bl0;
# Line 300 | Line 250
250            sprintf(hatext,"Bond constants missing for bond %d - %d of type %d - %d\n",bonds_ff.i12[i][0],bonds_ff.i12[i][1],
251             atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
252            fprintf(pcmlogfile,"%s\n",hatext);
303          //          fprintf(errfile,"%s\n",hatext);
253        }
254   L_10:
255        continue;
256            
257      }
258    }
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  }
259    if (Missing_constants == TRUE)
260    {
261        fprintf(pcmlogfile,"Error assigning constants in: %s\n",Struct_Title);
# Line 419 | Line 265
265            
266    return TRUE;    
267   }
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 }
268   /* ------------------------------------------  */
269   int find_bond(int ia, int ib)
270   {
# Line 450 | Line 283
283   }
284  
285   /* ------------------------------------------  */
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 /* ------------------------------------------  */
286   int is_delocalbond(int ia, int ib)
287   {
288      // test for delocalized bond
# Line 561 | Line 346
346          return( b < 0.0 ? -fabs(a) : fabs(a) );
347   }
348  
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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines