src/tools/dtools.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.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33 #include "smalloc.h"
  34 #include "strdb.h"
  35 #include "futil.h"
  36 #include "pdbio.h"
  37 #include "vec.h"
  38 #include "names.h"
  39 #include "disco.h"
  40
  41 char *edc_names[edcNR+1] = { "NONBND", "BOND", "DISRE", NULL };
  42
  43 bool newres(int i,t_atom atom[])
  44 {
  45   return ((i == 0) || (atom[i].resnr != atom[i-1].resnr));
  46 }
  47
  48 void pr_corr(FILE *log,t_correct *c)
  49 {
  50   fprintf(log,"Parameters for this disco run\n");
  51   fprintf(log,"maxnit        = %d\n",c->maxnit);
  52   fprintf(log,"nbcheck       = %d\n",c->nbcheck);
  53   fprintf(log,"nstprint      = %d\n",c->nstprint);
  54   fprintf(log,"nstranlist    = %d\n",c->nstranlist);
  55   fprintf(log,"ngrow         = %d\n",c->ngrow);
  56   fprintf(log,"bExplicit     = %s\n",BOOL(c->bExplicit));
  57   fprintf(log,"bChiral       = %s\n",BOOL(c->bChiral));
  58   fprintf(log,"bPep          = %s\n",BOOL(c->bPep));
  59   fprintf(log,"bDump         = %s\n",BOOL(c->bDump));
  60   fprintf(log,"bLowerOnly    = %s\n",BOOL(c->bLowerOnly));
  61   fprintf(log,"bRanlistFirst = %s\n",BOOL(c->bRanlistFirst));
  62   fprintf(log,"bCubic        = %s\n",BOOL(c->bCubic));
  63   fprintf(log,"bBox          = %s\n",BOOL(c->bBox));
  64   fprintf(log,"bCenter       = %s\n",BOOL(c->bCenter));
  65   fprintf(log,"ndist         = %d\n",c->ndist);
  66   fprintf(log,"npep          = %d\n",c->npep);
  67   fprintf(log,"nimp          = %d\n",c->nimp);
  68   fprintf(log,"lowdev        = %g\n",c->lodev);
  69   fflush(log);
  70 }
  71
  72 t_correct *init_corr(int maxnit,int nstprint,int nbcheck,int nstranlist,
  73                      int ngrow,bool bExplicit,bool bChiral,bool bPep,
  74                      bool bDump,real lowdev,bool bLowerOnly,
  75                      bool bRanlistFirst,bool bCubic,bool bBox,bool bCenter)
  76 {
  77   t_correct *c;
  78
  79   snew(c,1);
  80   c->maxnit        = maxnit;
  81   c->nstprint      = nstprint;
  82   c->nbcheck       = nbcheck;
  83   c->nstranlist    = nstranlist;
  84   c->bExplicit     = bExplicit;
  85   c->bChiral       = bChiral;
  86   c->bPep          = bPep;
  87   c->bDump         = bDump;
  88   c->lodev         = lowdev;
  89   c->maxdist       = 0;
  90   c->ndist         = 0;
  91   c->ngrow         = ngrow;
  92   c->bLowerOnly    = bLowerOnly;
  93   c->bRanlistFirst = bRanlistFirst;
  94   c->bCubic        = bCubic;
  95   c->bBox          = bBox;
  96   c->bCenter       = bCenter;
  97
  98   if (bRanlistFirst) {
  99     if (nstranlist != 0) {
 100       fprintf(stderr,
 101               "Warning: setting nstranlist to 0 when using ranlistfirst");
 102       c->nstranlist = 0;
 103     }
 104   }
 105   return c;
 106 }
 107
 108 void init_corr2(t_correct *c,int natom)
 109 {
 110   int i,n,ai;
 111
 112   if (c->ndist == 0)
 113     fatal_error(0,"Can not make tags without distances. %s, %d",
 114                 __FILE__,__LINE__);
 115
 116   /* Initiate the ip index */
 117   snew(c->ip,c->ndist);
 118
 119   /* Init boolean array */
 120   snew(c->bViol,natom);
 121
 122   /* Initiate omega array */
 123   snew(c->omega,c->npep);
 124
 125   /* Initiate improper array */
 126   srenew(c->idih,c->nimp+1);
 127
 128   /* Now make the tags */
 129   snew(c->tag,natom+1);
 130   ai = c->d[0].ai;
 131   n  = 0;
 132   for(i=0; (i<c->ndist); i++) {
 133     c->ip[i] = i;
 134     if (c->d[i].ai != ai) {
 135       /* New tag starts here, but skip over possible missing atoms */
 136       while ((n < natom) && (n<c->d[i].ai)) {
 137         n++;
 138         c->tag[n] = i;
 139         ai = c->d[i].ai;
 140       }
 141       /*     else
 142              fatal_error(0,"Too many atoms, or distances not sorted");*/
 143     }
 144   }
 145   if (n < natom) {
 146     while (n < natom) {
 147       n++;
 148       c->tag[n] = i;
 149     }
 150   }
 151   else
 152     fatal_error(0,"Too many atoms, or distances not sorted");
 153   if (debug)
 154     fprintf(debug,"There are %d tags for %d atoms\n",n,natom);
 155
 156   if (debug)
 157     for(i=0; (i<=natom); i++)
 158       fprintf(debug,"tag[%5d] = %d\n",i,c->tag[i]);
 159 }
 160
 161 void center_in_box(int natom,rvec xin[],matrix box,rvec xout[])
 162 {
 163   int  i,m;
 164   rvec xcm,dx;
 165
 166   /* Dump the optimization trajectory to an xtc file  */
 167   /* Put the whole thing in the center of the box 1st */
 168   clear_rvec(xcm);
 169   for(i=0; (i<natom); i++) {
 170     copy_rvec(xin[i],xout[i]);
 171     rvec_inc(xcm,xin[i]);
 172   }
 173   for(m=0; (m<DIM); m++)
 174     dx[m] = 0.5*box[m][m]-xcm[m]/natom;
 175   for(i=0; (i<natom); i++)
 176     rvec_inc(xout[i],dx);
 177 }
 178
 179 void define_peptide_bonds(FILE *log,t_atoms *atoms,t_correct *c)
 180 {
 181   int    i,npep,naa,nlist;
 182   char   **aa;
 183   t_dlist *dlist;
 184
 185   naa   = get_strings("aminoacids.dat",&aa);
 186   dlist = mk_dlist(log,atoms,&nlist,TRUE,TRUE,FALSE,0,1,naa,aa);
 187   for(i=0; (i<naa); i++)
 188     sfree(aa[i]);
 189   sfree(aa);
 190
 191   npep  = 0;
 192   snew(c->pepbond,nlist);
 193   for(i=0; (i<nlist); i++) {
 194     if (has_dihedral(edOmega,&dlist[i])) {
 195       c->pepbond[npep].ai = dlist[i].atm.minC;
 196       c->pepbond[npep].aj = dlist[i].atm.minO;
 197       c->pepbond[npep].ak = dlist[i].atm.N;
 198       c->pepbond[npep].al = dlist[i].atm.H;
 199       c->pepbond[npep].am = dlist[i].atm.Cn[1];
 200       npep++;
 201     }
 202   }
 203   c->npep = npep;
 204   if (debug)
 205     pr_dlist(debug,nlist,dlist,1.0);
 206   sfree(dlist);
 207   fprintf(log,"There are %d peptide bonds\n",npep);
 208 }
 209
 210 static char *aname(t_atoms *atoms,int i)
 211 {
 212   static char buf[32];
 213
 214   sprintf(buf,"%4s%3d:%4s",
 215           *(atoms->resname[atoms->atom[i].resnr]),
 216           atoms->atom[i].resnr+1,
 217           *(atoms->atomname[i]));
 218
 219   return buf;
 220 }
 221
 222 int find_atom(t_atoms *atoms,int resnr,int j0,char *nm)
 223 {
 224   int j,aa=-1;
 225
 226   for(j=j0; (j < atoms->nr) && (atoms->atom[j].resnr == resnr); j++)
 227     if (strcmp(*atoms->atomname[j],nm) == 0) {
 228       aa = j;
 229       break;
 230     }
 231   return aa;
 232 }
 233
 234 void define_impropers(FILE *log,t_atoms *atoms,t_correct *c)
 235 {
 236   typedef struct {
 237     char *res,*aa[4];
 238   } t_impdef;
 239
 240   t_impdef id[] = {
 241     { NULL,   { "CA", "N",  "C",   "CB"  } },
 242     { NULL,   { "N",  "CA", "H1",  "H3"  } },
 243     { "LYSH", { "NZ", "CE", "HZ2", "HZ1" } },
 244     { "LYSH", { "NZ", "CE", "HZ1", "HZ3" } },
 245     { "LEU",  { "CG", "CD2", "CD1", "CB" } },
 246     { "VAL",  { "CB", "CG2", "CG1", "CA" } },
 247     { "ILE",  { "CB", "CG1", "CG2", "CA" } },
 248     { "THR",  { "CB", "OG1", "CG2", "CA" } }
 249   };
 250 #define NID asize(id)
 251   int i,j,k,l,aa[4],nimp;
 252
 253   nimp = 0;
 254   for(i=0; (i<atoms->nres); i++) {
 255     /* Skip until residue number */
 256     for(j=0; (j<atoms->nr) && (atoms->atom[j].resnr != i); j++)
 257       ;
 258     for(k=0; (k<NID); k++) {
 259       if ((id[k].res == NULL) ||
 260           (strcmp(id[k].res,*atoms->resname[i]) == 0)) {
 261         /* This (i) is the right residue to look for this improper (k) */
 262         for(l=0; (l<4); l++)
 263           aa[l] = find_atom(atoms,i,j,id[k].aa[l]);
 264         if ((aa[0] != -1) && (aa[1] != -1) && (aa[2] != -1) && (aa[3] != -1)) {
 265           srenew(c->imp,nimp+1);
 266           c->imp[nimp].ai = aa[0];
 267           c->imp[nimp].aj = aa[1];
 268           c->imp[nimp].ak = aa[2];
 269           c->imp[nimp].al = aa[3];
 270           nimp++;
 271         }
 272       }
 273     }
 274   }
 275   c->nimp = nimp;
 276
 277   fprintf(log,"There are %d impropers\n",c->nimp);
 278   if (debug) {
 279     fprintf(debug,"Overview of improper dihedrals\n");
 280     for(i=0; (i<c->nimp); i++) {
 281       fprintf(debug,"  %s",aname(atoms,c->imp[i].ai));
 282       fprintf(debug,"  %s",aname(atoms,c->imp[i].aj));
 283       fprintf(debug,"  %s",aname(atoms,c->imp[i].ak));
 284       fprintf(debug,"  %s\n",aname(atoms,c->imp[i].al));
 285     }
 286   }
 287 }
 288
 289 void pr_dist(FILE *fp,bool bHeader,t_correct *c,int i)
 290 {
 291   real ideal=0;
 292
 293   if (bHeader)
 294     fprintf(fp,"#%4s%5s%10s%10s%10s\n","ai","aj","ideal","lb","ub");
 295   switch (c->d[i].cons_type) {
 296   case edcBOND:
 297     ideal = 5*(c->d[i].lb+c->d[i].ub);
 298     break;
 299   case edcNONBOND:
 300     ideal = 0;
 301     break;
 302   case edcDISRE:
 303     ideal = -1;
 304     break;
 305   default:
 306     fatal_error(0,"cons_type for distance %d = %d\n",i,c->d[i].cons_type);
 307   }
 308   fprintf(fp,"%5d%5d%10.5f%10.5f%10.5f\n",1+c->d[i].ai,1+c->d[i].aj,
 309           ideal,10*c->d[i].lb,10*c->d[i].ub);
 310 }
 311
 312 void pr_distances(FILE *fp,t_correct *c)
 313 {
 314   int i;
 315
 316   for(i=0; (i<c->ndist); i++)
 317     pr_dist(fp,(i==0),c,i);
 318 }
 319
 320 void add_dist(t_correct *c,int ai,int aj,real lb,real ideal,real ub,real w[])
 321 {
 322   int n = c->ndist;
 323
 324   if ((w[ai] != 0) || (w[aj] != 0)) {
 325     if (n == c->maxdist) {
 326       c->maxdist += 100;
 327       srenew(c->d,c->maxdist);
 328     }
 329     c->d[n].ai = ai;
 330     c->d[n].aj = aj;
 331     if (ideal > 0)
 332       c->d[n].cons_type = edcBOND;
 333     else if (ideal == 0.0)
 334       c->d[n].cons_type = edcNONBOND;
 335     else
 336       c->d[n].cons_type = edcDISRE;
 337     c->d[n].lb = lb;
 338     c->d[n].ub = ub;
 339     c->d[n].wi = w[ai]/(w[ai]+w[aj]);
 340     c->ndist++;
 341   }
 342 }
 343
 344 void read_dist(FILE *log,char *fn,int natom,t_correct *c)
 345 {
 346   FILE   *fp;
 347   char   buf[256];
 348   int    ai,aj,i,nline=0;
 349   double ideal,lb,ub;
 350   int    nnn[edcNR];
 351
 352   for(i=0; (i<edcNR); i++)
 353     nnn[i] = 0;
 354
 355   fp = ffopen(fn,"r");
 356   while (fgets(buf,255,fp)) {
 357     nline ++;
 358     if (buf[0] != '#') {
 359       if (sscanf(buf,"%d%d%lf%lf%lf",&ai,&aj,&ideal,&lb,&ub) != 5)
 360         fprintf(stderr,"Error in %s, line %d\n",fn,nline);
 361       else {
 362         add_dist(c,ai-1,aj-1,lb*0.1,ideal*0.1,ub*0.1,c->weight);
 363         nnn[c->d[c->ndist-1].cons_type]++;
 364       }
 365     }
 366   }
 367   ffclose(fp);
 368
 369   fprintf(log,"Read distances from file %s\n",fn);
 370   for(i=0; (i<edcNR); i++)
 371     fprintf(log,"There were %d %s distances\n",
 372             nnn[i],edc_names[i]);
 373 }
 374
 375 real *read_weights(char *fn,int natom)
 376 {
 377   t_atoms newatoms;
 378   int     i;
 379   rvec    *x;
 380   matrix  box;
 381   char    title[256];
 382   real    *w;
 383
 384   /* Read the weights from the occupancy field in the pdb file */
 385   snew(x,natom);
 386   init_t_atoms(&newatoms,natom,TRUE);
 387   read_pdb_conf(fn,title,&newatoms,x,box,FALSE);
 388   snew(w,newatoms.nr);
 389   for(i=0; (i<newatoms.nr); i++)
 390     w[i] = newatoms.pdbinfo[i].occup;
 391   free_t_atoms(&newatoms);
 392   sfree(x);
 393
 394   return w;
 395 }
 396
 397 void check_dist(FILE *log,t_correct *c)
 398 {
 399   int  i;
 400   real tol=0.001;
 401
 402   fprintf(log,"Checking distances for internal consistency\n");
 403   for(i=0; (i<c->ndist); i++) {
 404     if ((c->d[i].ub != 0) && ((c->d[i].ub - c->d[i].lb) < tol)) {
 405       pr_dist(log,TRUE,c,i);
 406     }
 407   }
 408 }
 409