src/tools/gendr.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  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_gendr_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include <string.h>
  34 #include "string2.h"
  35 #include "strdb.h"
  36 #include "typedefs.h"
  37 #include "macros.h"
  38 #include "copyrite.h"
  39 #include "smalloc.h"
  40 #include "statutil.h"
  41 #include "confio.h"
  42 #include "genhydro.h"
  43
  44 typedef struct {
  45   char *key;
  46   int  nexp;
  47   char **exp;
  48 } t_expansion;
  49
  50 t_expansion *read_expansion_map(char *fn,int *nexpand)
  51 {
  52   char        ibuf[12],buf[12][10];
  53   char        **ptr;
  54   t_expansion *exp;
  55   int         i,k,nexp,nn;
  56
  57   nexp=get_lines(fn,&ptr);
  58
  59   snew(exp,nexp);
  60   for(i=0; (i<nexp); i++) {
  61     /* Let scanf do the counting... */
  62     nn=sscanf(ptr[i],"%s%s%s%s%s%s%s%s%s%s%s",
  63               ibuf,buf[0],buf[1],buf[2],buf[3],buf[4],
  64               buf[5],buf[6],buf[7],buf[8],buf[9]);
  65     if (nn <= 1)
  66       break;
  67     exp[i].key=strdup(ibuf);
  68     exp[i].nexp=nn-1;
  69     snew(exp[i].exp,nn-1);
  70     for(k=0; (k<nn-1); k++)
  71       exp[i].exp[k]=strdup(buf[k]);
  72   }
  73   fprintf(stderr,"I found %d expansion mapping entries!\n",i);
  74
  75   /* Clean up */
  76   for(i=0; (i<nexp); i++)
  77     sfree(ptr[i]);
  78   sfree(ptr);
  79
  80   *nexpand=nexp;
  81   return exp;
  82 }
  83
  84 char **get_exp(int NEXP,t_expansion expansion[],char **ai,int *nexp)
  85 {
  86   int  i;
  87
  88   for(i=0; (i<NEXP); i++)
  89     if (strcmp(*ai,expansion[i].key) == 0) {
  90       *nexp=expansion[i].nexp;
  91       return expansion[i].exp;
  92     }
  93   *nexp=1;
  94
  95   return ai;
  96 }
  97
  98 int find_atom(char *ai,char *ri,
  99               int resi,int r0,
 100               int natoms,char ***aname,t_atom atom[],
 101               int linec)
 102 {
 103   int i;
 104
 105   /* Locate residue */
 106   for(i=0; (i<natoms) && (atom[i].resnr != resi); i++)
 107     ;
 108   if (i == natoms)
 109     return -1;
 110
 111   /* Compare atom names */
 112   for(   ; (i<natoms) && (atom[i].resnr == resi); i++)
 113     if (strcmp(*(aname[i]),ai) == 0)
 114       return i;
 115
 116   /* Not found?! */
 117   fprintf(stderr,"Warning: atom %s not found in res %s%d (line %d)\n",
 118           ai,ri,resi+r0,linec);
 119
 120   return -1;
 121 }
 122
 123 void conv_dr(FILE *in,FILE *out,char *map,t_atoms *atoms,int r0)
 124 {
 125   static char *format="%s%d%s%s%d%s%lf%lf";
 126   int    i,j,nc,nindex,ni,nj,nunres;
 127   int    atomi,atomj,resi,resj;
 128   char   **aiexp,**ajexp;
 129   char   *ai,*aj;
 130   char   ri[10],rj[10];
 131   char   buf[1024];
 132   double ub,lb;
 133   int    linec;
 134   int    NEXP;
 135   t_expansion *exp;
 136
 137   exp=read_expansion_map(map,&NEXP);
 138
 139   nc=0;
 140   nindex=0;
 141   nunres=0;
 142   snew(ai,10);
 143   snew(aj,10);
 144   fprintf(out,"[ distance_restraints ]\n");
 145   linec=1;
 146
 147   while (fgets2(buf,1023,in) != NULL) {
 148     /* Parse the input string. If your file format is different,
 149      * modify it here...
 150      * If your file contains no spaces but colon (:) for separators
 151      * it may be easier to just replace those by a space using a
 152      * text editor.
 153      */
 154     sscanf(buf,format,ri,&resi,ai,rj,&resj,aj,&lb,&ub);
 155     aiexp=get_exp(NEXP,exp,&ai,&ni);
 156     ajexp=get_exp(NEXP,exp,&aj,&nj);
 157
 158     /* Turn bounds into nm */
 159     ub*=0.1;
 160     lb*=0.1;
 161
 162     /* Subtract starting residue to match topology */
 163     resi-=r0;
 164     resj-=r0;
 165
 166     /* Test whether residue names match
 167      * Again, if there is a mismatch between GROMACS names
 168      * and your file (eg. HIS vs. HISH) it may be easiest to
 169      * use your text editor...
 170      */
 171     if (strcmp(*atoms->resname[resi],ri) != 0) {
 172       fprintf(stderr,"Warning resname in disres file %s%d, in tpx file %s%d\n",
 173               ri,resi+r0,*atoms->resname[resi],resi+r0);
 174       nunres++;
 175     }
 176     /* Residue j */
 177     else if (strcmp(*atoms->resname[resj],rj) != 0) {
 178       fprintf(stderr,"Warning resname in disres file %s%d, in tpx file %s%d\n",
 179               rj,resj+r0,*atoms->resname[resj],resj+r0);
 180       nunres++;
 181     }
 182     else {
 183       /* Here, both residue names match ! */
 184       for(i=0; (i<ni); i++) {
 185         if ((atomi=find_atom(aiexp[i],ri,resi,r0,atoms->nr,
 186                              atoms->atomname,atoms->atom,linec)) == -1)
 187           nunres++;
 188         else {
 189           /* First atom is found now... */
 190           for(j=0; (j<nj); j++) {
 191             if ((atomj=find_atom(ajexp[j],ri,resj,r0,atoms->nr,
 192                                  atoms->atomname,atoms->atom,linec)) == -1)
 193               nunres++;
 194             else {
 195               /* BOTH atoms are found now! */
 196               fprintf(out,"%5d  %5d  1  %5d  1  %8.3f  %8.3f  %8.3f  %8.3f\n",
 197                       1+atomi,1+atomj,nindex,ub,ub+0.1,0.0,0.0);
 198               nc++;
 199             }
 200           }
 201         }
 202       }
 203     }
 204     nindex++;
 205     linec++;
 206   }
 207   fprintf(stderr,"Total number of NOES: %d\n",nindex);
 208   fprintf(stderr,"Total number of restraints: %d\n",nc);
 209   fprintf(stderr,"Total number of unresolved atoms: %d\n",nunres);
 210   if (nunres+nc != nindex)
 211     fprintf(stderr,"Holy Cow! some lines have disappeared.\n");
 212 }
 213
 214 int main (int argc,char *argv[])
 215 {
 216   static char *desc[] = {
 217     "gendr generates a distance restraint entry for a gromacs topology",
 218     "from another format. The format of the input file must be:[BR]",
 219     "resnr-i resname-i atomnm-i resnr-j resname-j atomnm-j lower upper[BR]"  ,
 220     "where lower and upper are the distance bounds.",
 221     "The entries must be separated by spaces, but may be otherwise in",
 222     "free format. Some expansion of templates like MB -> HB1, HB2 is done",
 223     "but this is not really well tested."
 224   };
 225
 226   static int r0=1;
 227   t_pargs pa[] = {
 228     { "-r", FALSE, etINT, {&r0}, "starting residue number" }
 229   };
 230
 231   FILE        *in,*out;
 232   t_topology  *top;
 233
 234   t_filenm fnm[] = {
 235     { efTPX, "-s", NULL, ffREAD  },
 236     { efDAT, "-d", NULL, ffREAD  },
 237     { efITP, "-o", NULL, ffWRITE },
 238     { efDAT, "-m", "expmap", ffREAD }
 239   };
 240 #define NFILE asize(fnm)
 241
 242   CopyRight(stderr,argv[0]);
 243   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,asize(desc),desc,
 244                     0,NULL);
 245
 246   fprintf(stderr,"Will subtract %d from res numbers in %s\n",
 247           r0,ftp2fn(efDAT,NFILE,fnm));
 248
 249   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 250
 251   in=opt2FILE("-d",NFILE,fnm,"r");
 252   out=ftp2FILE(efITP,NFILE,fnm,"w");
 253   conv_dr(in,out,opt2fn("-m",NFILE,fnm),&(top->atoms),r0);
 254   fclose(in);
 255   fclose(out);
 256
 257   thanx(stdout);
 258
 259   return 0;
 260 }
 261
 262