src/tools/sas2mat.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 <string.h>
  40
  41 #include "sysstuff.h"
  42 #include "matio.h"
  43 #include "copyrite.h"
  44 #include "macros.h"
  45 #include "statutil.h"
  46 #include "smalloc.h"
  47
  48 int main(int argc,char *argv[])
  49 {
  50   const char *desc[] = {
  51     "sas2mat converts matrix data in [IT]raw[it] format to X PixMap format,",
  52     "which can be digested by xpm2ps to make nice plots.",
  53     "These [IT]raw[it] data may be generated by g_rms, do_dssp or your",
  54     "own program.[PAR]",
  55     "The program prompts the user for some parameters:[PAR]",
  56     "[TT]Enter nres, res0, nframes, dt, t0, nlevels:[tt][PAR]",
  57     "In this context nres is the number of residues, res0 the starting residue",
  58     "dt is the time step, t0 is the starting time, nlevels is the number",
  59     "of levels for coloring. By default a greyscale colormap is generated."
  60   };
  61   static gmx_bool   bCol=FALSE;
  62   static char   *title="Area (nm^2)";
  63   static real   ssmin=-1,ssmax=-1,t0=0,dt=1;
  64   static int    nres=1,nframes=1,r0=0,nlevels=20,nskip=0;
  65   t_pargs pa[] = {
  66     { "-col",     FALSE,  etBOOL, &bCol,
  67       "The user is prompted for rgb lower and upper values" },
  68     { "-min",     FALSE,  etREAL, &ssmin,
  69       "Lower values for the data, calculated from the data by default" },
  70     { "-max",     FALSE,  etREAL, &ssmax,
  71       "Upper values for the data, see above" },
  72     { "-title",   FALSE,  etSTR,  &title,
  73       "Title for the graph" },
  74     { "-nlevel",  FALSE,  etINT,  &nlevels,
  75       "Number of levels in graph" },
  76     { "-nres",    FALSE,  etINT,  &nres,
  77       "Number of residues (Y-axis)" },
  78     { "-nframes", FALSE,  etINT,  &nframes,
  79       "Number of frames (Y-axis)" },
  80     { "-res0",    FALSE,  etINT,  &r0,
  81       "Number of first residue" },
  82     { "-nskip",   FALSE,  etINT,  &nskip,
  83       "Number of frames to skip after every frame" },
  84     { "-dt",      FALSE,  etREAL, &dt,
  85       "Time between time frames" },
  86     { "-t0",      FALSE,  etREAL, &t0,
  87       "Time of first time frame" }
  88   };
  89
  90   FILE   *in,*out;
  91   int    i,j,k,ihi;
  92   double s;
  93   real   **ss,lo,hi,s1min,s1max;
  94   real   *resnr,*t;
  95   gmx_bool   bCheck=TRUE;
  96   t_rgb  rlo,rhi;
  97   t_filenm fnm[] = {
  98     { efOUT, "-f", "area", ffREAD },
  99     { efXPM, "-o", "sas",  ffWRITE }
 100   };
 101 #define NFILE asize(fnm)
 102
 103   /* If we want to read all frames nskip must be greater than zero */
 104   nskip += 1;
 105
 106   CopyRight(stderr,argv[0]);
 107
 108   parse_common_args(&argc,argv,PCA_BE_NICE,NFILE,fnm,asize(pa),pa,asize(desc),desc,
 109                     0,NULL);
 110
 111   snew(ss,nres);
 112   snew(resnr,nres);
 113   snew(t,nframes);
 114   for(i=0; (i<nframes); i++)
 115     t[i]=t0+i*dt;
 116   for(i=0; (i<nres); i++) {
 117     snew(ss[i],nframes);
 118   }
 119   in=ftp2FILE(efOUT,NFILE,fnm,"r");
 120   for(i=k=0; (i<nframes); i++) {
 121     for(j=0; (j<nres); j++) {
 122       fscanf(in,"%lf",&s);
 123       ss[j][k]=s;
 124     }
 125     if (!nskip || ((i % nskip) == 0))
 126       k++;
 127   }
 128   ffclose(in);
 129   nframes=k;
 130
 131   lo=10000;
 132   hi=0;
 133   for(j=0; (j<nres); j++) {
 134     /* Find lowest SAS value and subtract that from all occurrences */
 135     s1min=10000;
 136     s1max=0;
 137     for(i=0; (i<nframes); i++) {
 138       s1min=min(s1min,ss[j][i]);
 139       s1max=max(s1max,ss[j][i]);
 140     }
 141     printf("res %d: ssmin=%g, ssmax=%g, diff=%g\n",j,s1min,s1max,s1max-s1min);
 142     hi=max(hi,s1max);
 143     lo=min(lo,s1min);
 144   }
 145   printf("Lowest and Highest SAS value: %g %g\n",lo,hi);
 146
 147   if (ssmin == -1)
 148     ssmin=lo;
 149   if (ssmax == -1)
 150     ssmax=hi;
 151
 152   /*
 153     hi=ssmax-ssmin;
 154     for(j=0; (j<nres); j++) {
 155     for(i=0; (i<nframes); i++)
 156     ss[j][i]-=ssmin;
 157     }
 158     */
 159
 160   /* ihi=hi; */
 161   rhi.r=0,rhi.g=0,rhi.b=0;
 162   rlo.r=1,rlo.g=1,rlo.b=1;
 163   if (bCol) {
 164     printf("Color entries:\n""drlo glo blo rhi ghi bhi\n");
 165     scanf("%f%f%f%f%f%f",&rlo.r,&rlo.g,&rlo.b,&rhi.r,&rhi.g,&rhi.b);
 166   }
 167   /*
 168   write_mapfile(ftp2fn(efMAP,NFILE,fnm),&nlevels,rlo,rhi,ssmin,ssmax);
 169   */
 170
 171   for(i=0;i<nres;i++)
 172     resnr[i]=i+1;
 173   out=ftp2FILE(efXPM,NFILE,fnm,"w");
 174   /*
 175   write_matrix(out,nres,nframes,resnr,t,ss,NULL,title,0,hi,nlevels);
 176   */
 177   write_xpm(out,0,"????","????","Time (ps)","Residue",
 178             nres,nframes,resnr,t,ss,ssmin,ssmax,rlo,rhi,&nlevels);
 179
 180   thanx(stderr);
 181
 182   return 0;
 183 }