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