src/tools/gendr.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <math.h>
  41 #include <string.h>
  42 #include "string2.h"
  43 #include "strdb.h"
  44 #include "typedefs.h"
  45 #include "macros.h"
  46 #include "copyrite.h"
  47 #include "smalloc.h"
  48 #include "statutil.h"
  49 #include "confio.h"
  50 #include "calch.h"
  51
  52 typedef struct {
  53   char *key;
  54   int  nexp;
  55   char **exp;
  56 } t_expansion;
  57
  58 t_expansion *read_expansion_map(char *fn,int *nexpand)
  59 {
  60   char        ibuf[12],buf[12][10];
  61   char        **ptr;
  62   t_expansion *exp;
  63   int         i,k,nexp,nn;
  64
  65   nexp=get_file(fn,&ptr);
  66
  67   snew(exp,nexp);
  68   for(i=0; (i<nexp); i++) {
  69     /* Let scanf do the counting... */
  70     nn=sscanf(ptr[i],"%s%s%s%s%s%s%s%s%s%s%s",
  71               ibuf,buf[0],buf[1],buf[2],buf[3],buf[4],
  72               buf[5],buf[6],buf[7],buf[8],buf[9]);
  73     if (nn <= 1)
  74       break;
  75     exp[i].key=strdup(ibuf);
  76     exp[i].nexp=nn-1;
  77     snew(exp[i].exp,nn-1);
  78     for(k=0; (k<nn-1); k++)
  79       exp[i].exp[k]=strdup(buf[k]);
  80   }
  81   fprintf(stderr,"I found %d expansion mapping entries!\n",i);
  82
  83   /* Clean up */
  84   for(i=0; (i<nexp); i++)
  85     sfree(ptr[i]);
  86   sfree(ptr);
  87
  88   *nexpand=nexp;
  89   return exp;
  90 }
  91
  92 char **get_exp(int NEXP,t_expansion expansion[],char **ai,int *nexp)
  93 {
  94   int  i;
  95
  96   for(i=0; (i<NEXP); i++)
  97     if (strcmp(*ai,expansion[i].key) == 0) {
  98       *nexp=expansion[i].nexp;
  99       return expansion[i].exp;
 100     }
 101   *nexp=1;
 102
 103   return ai;
 104 }
 105
 106 int find_atom(char *ai,char *ri,
 107               int resi,int r0,
 108               int natoms,char ***aname,t_atom atom[],
 109               int linec,bool bVerbose)
 110 {
 111   int i;
 112
 113   /* Locate residue */
 114   for(i=0; (i<natoms) && (atom[i].resnr != resi); i++)
 115     ;
 116   if (i == natoms)
 117     return -1;
 118
 119   /* Compare atom names */
 120   for(   ; (i<natoms) && (atom[i].resnr == resi); i++)
 121     if (strcmp(*(aname[i]),ai) == 0)
 122       return i;
 123
 124   /* Not found?! */
 125   if (bVerbose)
 126     fprintf(stderr,"Warning: atom %s not found in res %s%d (line %d)\n",
 127             ai,ri ? ri : "",resi+r0,linec);
 128
 129   return -1;
 130 }
 131
 132 void conv_dr(FILE *in,FILE *out,char *map,t_atoms *atoms,int r0,bool bXplor,
 133              bool bVerbose)
 134 {
 135   static char *format="%s%d%s%s%d%s%lf%lf";
 136   static char *xplorformat="%d%s%d%s";
 137   bool   bOK;
 138   int    i,j,nc,nindex,ni,nj,nunres;
 139   int    atomi,atomj,resi,resj;
 140   char   **aiexp,**ajexp;
 141   char   *ai,*aj;
 142   char   *ri,*rj;
 143   char   buf[1024];
 144   double ub,lb;
 145   int    linec;
 146   int    NEXP;
 147   t_expansion *exp;
 148
 149   exp=read_expansion_map(map,&NEXP);
 150
 151   nc=0;
 152   nindex=0;
 153   nunres=0;
 154   snew(ai,10);
 155   snew(aj,10);
 156   fprintf(out,"[ distance_restraints ]\n");
 157   linec=1;
 158
 159   if (bXplor) {
 160     ri = rj = NULL;
 161   }
 162   else {
 163     snew(ri,16);
 164     snew(rj,16);
 165   }
 166   while (fgets2(buf,1023,in) != NULL) {
 167     /* Parse the input string. If your file format is different,
 168      * modify it here...
 169      * If your file contains no spaces but colon (:) for separators
 170      * it may be easier to just replace those by a space using a
 171      * text editor.
 172      */
 173     if (bXplor) {
 174       bOK = (sscanf(buf,xplorformat,&resi,ai,&resj,aj) == 4);
 175       /* Cut atomnames at 4 characters */
 176       if (strlen(ai) >= 4)
 177         ai[4] = '\0';
 178       if (strlen(aj) >= 4)
 179         aj[4] = '\0';
 180       ub = 5.0;
 181       lb = 2.0;
 182     }
 183     else {
 184       bOK = (sscanf(buf,format,ri,&resi,ai,rj,&resj,aj,&lb,&ub) == 8);
 185     }
 186     if (bOK) {
 187       aiexp=get_exp(NEXP,exp,&ai,&ni);
 188       ajexp=get_exp(NEXP,exp,&aj,&nj);
 189
 190       /* Turn bounds into nm */
 191       ub*=0.1;
 192       lb*=0.1;
 193
 194       /* Subtract starting residue to match topology */
 195       resi-=r0;
 196       resj-=r0;
 197
 198       /* Test whether residue names match
 199        * Again, if there is a mismatch between GROMACS names
 200        * and your file (eg. HIS vs. HISH) it may be easiest to
 201        * use your text editor...
 202        */
 203
 204       if (!bXplor) {
 205         bOK = (strcmp(*atoms->resname[resi],ri) == 0);
 206         if (!bOK) {
 207           fprintf(stderr,"Warning resname in disres file %s%d, in tpx file %s%d\n",
 208                   ri,resi+r0,*atoms->resname[resi],resi+r0);
 209           nunres++;
 210         }
 211         else {
 212           /* Residue j */
 213           bOK = (strcmp(*atoms->resname[resj],rj) != 0);
 214           if (!bOK) {
 215             fprintf(stderr,"Warning resname in disres file %s%d, in tpx file %s%d\n",
 216                     rj,resj+r0,*atoms->resname[resj],resj+r0);
 217             nunres++;
 218           }
 219         }
 220       }
 221       if (bOK) {
 222         /* Here, both residue names match ! */
 223         for(i=0; (i<ni); i++) {
 224           if ((atomi=find_atom(aiexp[i],ri,resi,r0,atoms->nr,
 225                                atoms->atomname,atoms->atom,linec,bVerbose)) == -1)
 226             nunres++;
 227           else {
 228             /* First atom is found now... */
 229             for(j=0; (j<nj); j++) {
 230               if ((atomj=find_atom(ajexp[j],rj,resj,r0,atoms->nr,
 231                                    atoms->atomname,atoms->atom,linec,bVerbose)) == -1)
 232                 nunres++;
 233               else {
 234                 /* BOTH atoms are found now! */
 235                 fprintf(out,"%5d  %5d  1  %5d  1  %8.3f  %8.3f  %8.3f  %8.3f\n",
 236                         1+atomi,1+atomj,nindex,lb,ub,0.0,0.0);
 237                 nc++;
 238               }
 239             }
 240           }
 241         }
 242       }
 243       nindex++;
 244     }
 245     linec++;
 246   }
 247   fprintf(stderr,"Total number of NOES: %d\n",nindex);
 248   fprintf(stderr,"Total number of restraints: %d\n",nc);
 249   fprintf(stderr,"Total number of unresolved atoms: %d\n",nunres);
 250   if (nunres+nc != nindex)
 251     fprintf(stderr,"Holy Cow! some lines have disappeared.\n");
 252 }
 253
 254 int main (int argc,char *argv[])
 255 {
 256   const char *desc[] = {
 257     "gendr generates a distance restraint entry for a gromacs topology",
 258     "from another format. The format of the input file must be:[BR]",
 259     "resnr-i resname-i atomnm-i resnr-j resname-j atomnm-j lower upper[BR]"  ,
 260     "where lower and upper are the distance bounds.",
 261     "The entries must be separated by spaces, but may be otherwise in",
 262     "free format. Some expansion of templates like MB -> HB1, HB2 is done",
 263     "but this is not really well tested."
 264   };
 265   const char *bugs[] = {
 266     "This program is not well tested. Use at your own risk."
 267   };
 268
 269   static int  r0       = 1;
 270   static bool bXplor   = FALSE;
 271   static bool bVerbose = FALSE;
 272   t_pargs pa[] = {
 273     { "-r",     FALSE, etINT,  {&r0},       "starting residue number" },
 274     { "-xplor", FALSE, etBOOL, {&bXplor},   "Use xplor format for input" },
 275     { "-v",     FALSE, etBOOL, {&bVerbose}, "Be loud and noisy" }
 276   };
 277
 278   FILE        *in,*out;
 279   t_topology  *top;
 280
 281   t_filenm fnm[] = {
 282     { efTPX, "-s", NULL, ffREAD  },
 283     { efDAT, "-d", NULL, ffREAD  },
 284     { efITP, "-o", NULL, ffWRITE },
 285     { efDAT, "-m", "expmap", ffREAD }
 286   };
 287 #define NFILE asize(fnm)
 288
 289   CopyRight(stderr,argv[0]);
 290   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,asize(desc),desc,
 291                     asize(bugs),bugs);
 292
 293   fprintf(stderr,"******************* WARNING *****************************\n");
 294   fprintf(stderr,"*** Use at your own risk. When in doubt check the source.\n");
 295   fprintf(stderr,"*** Hang on: check the source anyway.\n");
 296   fprintf(stderr,"******************* WARNING *****************************\n");
 297
 298   fprintf(stderr,"Will subtract %d from res numbers in %s\n",
 299           r0,ftp2fn(efDAT,NFILE,fnm));
 300
 301   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 302
 303   in  = opt2FILE("-d",NFILE,fnm,"r");
 304   out = ftp2FILE(efITP,NFILE,fnm,"w");
 305   conv_dr(in,out,opt2fn("-m",NFILE,fnm),&(top->atoms),r0,bXplor,bVerbose);
 306   ffclose(in);
 307   ffclose(out);
 308
 309   thanx(stderr);
 310
 311   return 0;
 312 }
 313
 314