src/tools/bondlist.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_bondlist_c = "$Id$";
  30
  31 typedef struct list {
  32   char *residue;
  33   struct atombond {
  34     char *ai,*aj;
  35     real dist;
  36   } *atombond;
  37   struct atomangle {
  38     char *ai,*aj;
  39     real angle;
  40   } *atomangle;
  41   int natombond;
  42   int natomangle;
  43   struct list *next;
  44 } bondlist;
  45
  46 static bondlist *lst=NULL;
  47
  48 static void newresidue(char *residue)
  49 {
  50   bondlist *p;
  51
  52   snew(p,1);
  53   p->residue    = strdup(residue);
  54   p->atombond   = NULL;
  55   p->atomangle  = NULL;
  56   p->natombond  = 0;
  57   p->natomangle = 0;
  58   p->next       = lst;
  59   lst           = p;
  60   if (debug)
  61     fprintf(debug,"initialised element for residue %s\n",residue);
  62 }
  63
  64 void read_O_dist(void)
  65 {
  66 #define BLEN 255
  67   static   char *fn= "refi_aa.dat";
  68   FILE     *fp;
  69   char     buf[BLEN+1],buf_ai[9],buf_aj[9],junk1[15],junk2[15];
  70   int      nres = 0;
  71   real     dist,junk3;
  72
  73   fp = libopen(fn);
  74   fprintf(stderr,"Going to read %s\n",fn);
  75   fprintf(stderr,"Take note that the angle O C N+ will not be read from the");
  76   fprintf(stderr," correct residue\nThis is not a problem in the current");
  77   fprintf(stderr,"(980927) param-file while \nthis angle is const. 123.000\n");
  78   while (fgets2(buf,BLEN,fp) != NULL) {
  79     switch (buf[0]) {
  80     case '#':
  81     case '.':
  82     case 't':
  83       continue;
  84       break;
  85     case 'r':
  86       sscanf(buf,"%s%s",buf_ai,buf_aj);
  87       if (strcmp(buf_ai,"residue") == 0 ) {
  88         nres++;
  89         /* O only uses HIS, Gromos has HISB... */
  90         if (strcmp(buf_aj,"HIS") == 0) strcpy(buf_aj,"HISB");
  91         newresidue(buf_aj);
  92       }
  93       break;
  94     case 'b':
  95       if ( buf[5] == 'd' ) {
  96         sscanf(buf,"%*s%s%s%lf",buf_ai,buf_aj,&dist);
  97         /*fprintf(stderr,"1: %s %s %g\n",buf_ai,buf_aj,dist);*/
  98         if (strcmp(buf_ai,"C-") == 0) strcpy(buf_ai,"C");
  99         if (strcmp(buf_ai,"CA-") == 0) strcpy(buf_ai,"CA");
 100         if (strcmp(buf_ai,"N+") == 0) strcpy(buf_ai,"N");
 101         if (strcmp(buf_aj,"C-") == 0) strcpy(buf_aj,"C");
 102         if (strcmp(buf_aj,"CA-") == 0) strcpy(buf_aj,"CA");
 103         if (strcmp(buf_aj,"N+") == 0) strcpy(buf_aj,"N");
 104         /*fprintf(stderr,"2: %s %s %g\n",buf_ai,buf_aj,dist);*/
 105         lst->natombond++;
 106         lst->atombond = srenew(lst->atombond,lst->natombond);
 107         lst->atombond[lst->natombond-1].ai = strdup(buf_ai);
 108         lst->atombond[lst->natombond-1].aj = strdup(buf_aj);
 109         lst->atombond[lst->natombond-1].dist = (dist/10.0);
 110       }
 111       if ( buf[5] == 'a' ) {
 112         sscanf(buf,"%*s%s%*s%s%lf",buf_ai,buf_aj,&dist);
 113         /*fprintf(stderr,"1: %s %s %g\n",buf_ai,buf_aj,dist);*/
 114         if (strcmp(buf_ai,"C-") == 0) strcpy(buf_ai,"C");
 115         if (strcmp(buf_ai,"CA-") == 0) strcpy(buf_ai,"CA");
 116         if (strcmp(buf_ai,"N+") == 0) strcpy(buf_ai,"N");
 117         if (strcmp(buf_aj,"C-") == 0) strcpy(buf_aj,"C");
 118         if (strcmp(buf_aj,"CA-") == 0) strcpy(buf_aj,"CA");
 119         if (strcmp(buf_aj,"N+") == 0) strcpy(buf_aj,"N");
 120         /*fprintf(stderr,"2: %s %s %g\n",buf_ai,buf_aj,dist);*/
 121         lst->natomangle++;
 122         lst->atomangle = srenew(lst->atomangle,lst->natomangle);
 123         lst->atomangle[lst->natomangle-1].ai = strdup(buf_ai);
 124         lst->atomangle[lst->natomangle-1].aj = strdup(buf_aj);
 125         lst->atomangle[lst->natomangle-1].angle = dist;
 126       }
 127       break;
 128     default:
 129       fprintf(stderr,"Ooops! Unexpected first character in line:\n%s",buf);
 130       continue;
 131       break;
 132     }
 133   }
 134   fclose(fp);
 135   fprintf(stderr,"Read distances for %d residues from file %s \n",nres,fn);
 136 }
 137
 138 static real do_specialbond(char *name1,char *name2,char *res)
 139 {
 140   bondlist *p;
 141   int      i;
 142
 143   if (debug)
 144     fprintf(debug,"Looking for dist between %s %s in res. %s\n",
 145             name1,name2,res);
 146
 147   for ( p = lst; p != NULL; p = p->next)
 148     if (strcmp(res,p->residue) == 0)
 149       break;
 150
 151   /* Here either p == NULL or found */
 152   for (i=0; p!= NULL && (i<p->natombond) ; i++) {
 153     if (((strcmp(name1,p->atombond[i].ai) == 0) &&
 154          (strcmp(name2,p->atombond[i].aj) == 0)) ||
 155         ((strcmp(name1,p->atombond[i].aj) == 0) &&
 156          (strcmp(name2,p->atombond[i].ai) == 0))) {
 157       if (debug) {
 158         fprintf(debug,"Found dist. %s %s %s",name1,name2,res);
 159         fprintf(debug,"   %g\n",p->atombond[i].dist);
 160       }
 161       return p->atombond[i].dist;
 162     }
 163   }
 164   if (debug)
 165     fprintf(debug,"Didn't find dist. %s %s %s\n",name1,name2,res);
 166
 167   return 0.0;
 168 }
 169
 170 static real do_specialangle(char *name1,char *name2,char *res)
 171 {
 172
 173   bondlist *p;
 174   int      i;
 175
 176   if (debug)
 177     fprintf(debug,"Looking for angle between %s %s in res. %s\n",
 178             name1,name2, res);
 179
 180   for ( p = lst; p != NULL; p = p->next)
 181     if (strcmp(res,p->residue) == 0)
 182       break;
 183
 184   /* Here either p == NULL or found */
 185   for (i=0; p != NULL && (i<p->natomangle) ; i++) {
 186     if (((strcmp(name1,p->atomangle[i].ai) == 0) &&
 187          (strcmp(name2,p->atomangle[i].aj) == 0)) ||
 188         ((strcmp(name1,p->atomangle[i].aj) == 0) &&
 189          (strcmp(name2,p->atomangle[i].ai) == 0))) {
 190       if (debug) {
 191         fprintf(debug,"Found angle. %s %s %s",name1,name2,res);
 192         fprintf(debug,"   %g\n",p->atomangle[i].angle);
 193       }
 194       return p->atomangle[i].angle;
 195     }
 196   }
 197   if (debug)
 198     fprintf(debug,"Didn't find angle %s %s %s\n",name1,name2,res);
 199
 200   return 0.0;
 201 }
 202
 203 real lookup_bondlength(int ai,int aj,t_ilist ilist[],
 204                        t_iparams iparams[],bool bFail,t_atoms *atoms)
 205 {
 206   int  dodist[] = { F_BONDS, F_MORSE, F_SHAKE, F_G96BONDS };
 207   int  i,j,type,ftype,aai,aaj,nra;
 208   real blen=NOBOND;
 209   real spclen;
 210
 211   if (debug)
 212     fprintf(debug,"call do_spec with %s (%d) %s (%d) res %s \n",
 213             *atoms->atomname[ai],ai,*atoms->atomname[aj],aj,
 214             *(atoms->resname[atoms->atom[ai].resnr]));
 215
 216   spclen = do_specialbond(*(atoms->atomname[ai]),*(atoms->atomname[aj]),
 217                           *(atoms->resname[atoms->atom[aj].resnr]));
 218
 219   /* Comment; if you're using distances from e.g. O, it is non-trivial
 220      which residuename (ai or aj) you call do_special with, in order
 221      to get the peptide bonds correct. */
 222   if (spclen != 0) {
 223     if (debug)
 224       fprintf(debug,"%s %d %d %g\n",*atoms->resname[atoms->atom[ai].resnr],
 225               ai,aj,spclen);
 226     blen = spclen;
 227   }
 228   else {
 229     if (debug)
 230       fprintf(debug,"No spclen found for %s (%d) %s (%d) res %s \n",
 231               *atoms->atomname[ai],ai,*atoms->atomname[aj],aj,
 232               *(atoms->resname[atoms->atom[ai].resnr]));
 233     for(j=0; (j<asize(dodist)) && (blen == NOBOND); j++) {
 234       ftype = dodist[j];
 235       nra   = interaction_function[ftype].nratoms;
 236
 237       for(i=0; (i<ilist[ftype].nr) && (blen == NOBOND); i+=nra+1) {
 238         type  = ilist[ftype].iatoms[i];
 239         aai    = ilist[ftype].iatoms[i+1];
 240         aaj    = ilist[ftype].iatoms[i+2];
 241
 242         if (((aai == ai) && (aaj == aj)) || ((aaj == ai) && (aai == aj))) {
 243           switch (ftype) {
 244           case F_G96BONDS:
 245             blen = sqrt(iparams[type].harmonic.rA);
 246             break;
 247           case F_BONDS:
 248             blen = iparams[type].harmonic.rA;
 249             break;
 250           case F_MORSE:
 251             blen = iparams[type].morse.b0;
 252             break;
 253           case F_SHAKE:
 254             blen = iparams[type].shake.dA;
 255             break;
 256           default:
 257             break;
 258           }
 259         }
 260       }
 261     }
 262   }
 263   if (blen == NOBOND) {
 264     if (bFail)
 265       fatal_error(0,"No bond between atoms %d and %d\n",ai,aj);
 266     else
 267       return NOBOND;
 268   }
 269
 270   return blen*10;
 271   /* This should be the only conversion from nanometer to Angstrom */
 272 }
 273
 274 real lookup_angle(int ai,int aj,int ak,t_ilist ilist[],
 275                   t_iparams iparams[],t_atoms *atoms)
 276 {
 277   int  ft[] = { F_ANGLES, F_G96ANGLES };
 278   int  i,j,type,ff,ftype,aai,aaj,aak,nra;
 279   real angle;
 280
 281   /* First check the Engh & Huber database */
 282   angle = DEG2RAD*do_specialangle(*atoms->atomname[ai],*atoms->atomname[ak],
 283                                   *atoms->resname[atoms->atom[ak].resnr]);
 284   /* Now check the topology */
 285   for(ff=0; (ff<asize(ft)) && (angle == 0); ff++) {
 286     ftype = ft[ff];
 287     nra   = interaction_function[ftype].nratoms;
 288
 289     for(i=0; (i<ilist[ftype].nr) && (angle == 0); i+=nra+1) {
 290       type  = ilist[ftype].iatoms[i];
 291       aai   = ilist[ftype].iatoms[i+1];
 292       aaj   = ilist[ftype].iatoms[i+2];
 293       aak   = ilist[ftype].iatoms[i+3];
 294
 295       if (((aai == ai) && (aaj == aj) && (aak == ak)) ||
 296           ((aak == ai) && (aaj == aj) && (aai == ak))) {
 297         if (ftype == F_ANGLES)
 298           angle = DEG2RAD*iparams[type].harmonic.rA;
 299         else if (ftype == F_G96ANGLES)
 300           angle = acos(iparams[type].harmonic.rA);
 301         else
 302           fatal_error(0,"Unknown angletype %s in %s, line %d",
 303                       interaction_function[ftype].longname,__FILE__,__LINE__);
 304       }
 305     }
 306   }
 307   if (angle == 0) {
 308     fprintf(stderr,
 309             "Warning: no known angle between atoms %d, %d, %d. Using 109.47\n",
 310             ai,aj,ak);
 311     angle = 109.47*DEG2RAD;
 312   }
 313   return angle;
 314 }
 315
 316 real angle_length(int ai,int aj,int ak,real theta,
 317                   t_ilist ilist[],t_iparams iparams[],t_atoms *atoms)
 318 {
 319   real rij,rjk;
 320
 321   rij = lookup_bondlength(ai,aj,ilist,iparams,TRUE,atoms);
 322   rjk = lookup_bondlength(aj,ak,ilist,iparams,TRUE,atoms);
 323
 324   if (debug)
 325     fprintf(debug,"angle_length uses %g %g and angle %g\n",rij,rjk,theta);
 326
 327   return sqrt(rij*rij+rjk*rjk-2.0*rij*rjk*cos(DEG2RAD*theta));
 328 }
 329
 330 /*
 331   static void special_bonds(t_dist *d,t_atoms *atoms,
 332   char *res,int nbonds,t_bonddef bdef[])
 333   {
 334   int i,j,k,resnr,na,ai,aj,ndist;
 335
 336   resnr = -1;
 337   for(i=0; (i<atoms->nr); i++) {
 338     if (strcmp(*atoms->resname[atoms->atom[i].resnr],res) == 0) {
 339       resnr = atoms->atom[i].resnr;
 340       for(k=0; (k<nbonds); k++) {
 341         ai = aj = -1;
 342         for(j=i; ((j<atoms->nr) && (atoms->atom[j].resnr == resnr)); j++) {
 343           if (strcmp(bdef[k].ai,*atoms->atomname[j]) == 0)
 344             if (ai != -1)
 345               fatal_error(0,"Atom %s multiply defined in res %s %d",
 346                           bdef[k].ai,res,resnr);
 347             else
 348               ai = j;
 349           if (strcmp(bdef[k].aj,*atoms->atomname[j]) == 0)
 350             if (aj != -1)
 351               fatal_error(0,"Atom %s multiply defined in res %s %d",
 352                           bdef[k].aj,res,resnr);
 353             else
 354               aj = j;
 355         }
 356         if ((ai != -1) && (aj != -1)) {
 357           lb=(1.0-pep_margin)*bdef[k].len;
 358           ub=(1.0+pep_margin)*bdef[k].len;
 359           if (((weight[ai] != 0.0) || (weight[aj] != 0.0)) &&
 360               (!dist_set(d,natoms,ai,aj))) {
 361             set_dist(d,natoms,ai,aj,lb,ub,len);
 362             ndist++;
 363           }
 364         }
 365       }
 366       } */
 367     /* Optimizations ... */
 368 /*  }
 369    fprintf(stderr,"There were %d special distances defined (GROMOS-override)",
 370            ndist);
 371  }
 372 */
 373
 374
 375 void pdih_lengths(int ai,int aj,int ak,int al,
 376                   t_ilist ilist[],t_iparams iparams[],
 377                   real *lb,real *ub,t_atoms *atoms)
 378 {
 379   /* Returns the length corresponding to a cis dihedral */
 380   real rij,rjk,rkl,rik;
 381   real th1,th2,th3;
 382   real half_pi = M_PI*0.5;
 383
 384   rij = lookup_bondlength(ai,aj,ilist,iparams,TRUE,atoms);
 385   rjk = lookup_bondlength(aj,ak,ilist,iparams,TRUE,atoms);
 386   rkl = lookup_bondlength(ak,al,ilist,iparams,TRUE,atoms);
 387   th1 = lookup_angle(ai,aj,ak,ilist,iparams,atoms);
 388   th2 = lookup_angle(aj,ak,al,ilist,iparams,atoms);
 389
 390   /* Compute distance from i to k using law of cosines */
 391   rik = sqrt(rij*rij+rjk*rjk-2.0*rij*rjk*cos(th1));
 392
 393   /* Compute angle th3 using law of sines */
 394   th3 = asin(rij*sin(th1)/rik);
 395
 396   /* Compute cis length using law of cosines */
 397   *lb = sqrt(rik*rik+rkl*rkl-2.0*rik*rkl*cos(th2-th3));
 398
 399   /* Compute trans length using law of cosines */
 400   *ub = sqrt(rik*rik+rkl*rkl-2.0*rik*rkl*cos(th2+th3));
 401 }
 402
 403 real idih_lengths(int ai,int aj,int ak,int al,
 404                   t_ilist ilist[],t_iparams iparams[],t_atoms *atoms)
 405 {
 406   /* Returns the length corresponding to a cis dihedral */
 407   real rij,rjk,rkl,rik;
 408   real th1,th2,th3,cis;
 409   real half_pi = M_PI*0.5;
 410   real lb;
 411
 412   lb = 0.0;
 413
 414   if ((rij = lookup_bondlength(ai,aj,ilist,iparams,FALSE,atoms)) == NOBOND)
 415     return lb;
 416   if ((rjk = lookup_bondlength(aj,ak,ilist,iparams,FALSE,atoms)) == NOBOND)
 417     return lb;
 418   if ((rkl = lookup_bondlength(ak,al,ilist,iparams,FALSE,atoms)) == NOBOND)
 419     return lb;
 420
 421   if (debug)
 422     fprintf(debug,"Found an improper: %d %d %d %d\n",ai,aj,ak,al);
 423   th1 = lookup_angle(ai,aj,ak,ilist,iparams,atoms);
 424   th2 = lookup_angle(aj,ak,al,ilist,iparams,atoms);
 425
 426   /* Compute distance from i to k using law of cosines */
 427   rik = sqrt(rij*rij+rjk*rjk-2.0*rij*rjk*cos(th1));
 428
 429   /* Compute angle th3 using law of sines */
 430   th3 = asin(rij*sin(th1)/rik);
 431
 432   /* Compute cis length using law of cosines */
 433   return sqrt(rik*rik+rkl*rkl-2.0*rik*rkl*cos(th2-th3));
 434 }