src/contrib/g_anavel.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  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_com_c = "$Id$";
  30
  31 #include "sysstuff.h"
  32 #include "smalloc.h"
  33 #include "macros.h"
  34 #include "statutil.h"
  35 #include "random.h"
  36 #include "names.h"
  37 #include "matio.h"
  38 #include "physics.h"
  39 #include "vec.h"
  40 #include "futil.h"
  41 #include "copyrite.h"
  42 #include "xvgr.h"
  43 #include "string2.h"
  44 #include "rdgroup.h"
  45 #include "tpxio.h"
  46
  47 int main(int argc,char *argv[])
  48 {
  49   static char *desc[] = {
  50     "g_anavel computes temperature profiles in a sample. The sample",
  51     "can be analysed radial, i.e. the temperature as a function of",
  52     "distance from the center, cylindrical, i.e. as a function of distance",
  53     "from the vector (0,0,1) through the center of the box, or otherwise",
  54     "(will be specified later)"
  55   };
  56   t_filenm fnm[] = {
  57     { efTRN,  "-f",  NULL, ffREAD },
  58     { efTPX,  "-s",  NULL, ffREAD },
  59     { efXPM,  "-o", "xcm", ffWRITE }
  60   };
  61 #define NFILE asize(fnm)
  62
  63   static int  mode = 0,   nlevels = 10;
  64   static real tmax = 300, xmax    = -1;
  65   t_pargs pa[] = {
  66     { "-mode",    FALSE, etINT,  {&mode},    "mode" },
  67     { "-nlevels", FALSE, etINT,  {&nlevels}, "number of levels" },
  68     { "-tmax",    FALSE, etREAL, {&tmax},    "max temperature in output" },
  69     { "-xmax",    FALSE, etREAL, {&xmax},    "max distance from center" }
  70   };
  71
  72   FILE       *fp;
  73   int        *npts,nmax;
  74   int        status;
  75   int        i,j,idum,step,nframe=0,natoms,index;
  76   real       t,temp,rdum,hboxx,hboxy,scale,xnorm;
  77   real       **profile=NULL;
  78   real       *t_x=NULL,*t_y,hi=0;
  79   rvec       *x,*v;
  80   t_topology *top;
  81   int        d,m,n;
  82   matrix     box;
  83   atom_id    *sysindex;
  84   bool       bHaveV,bReadV;
  85   t_rgb      rgblo = { 0, 0, 1 },rgbhi = { 1, 0, 0 };
  86
  87
  88   CopyRight(stderr,argv[0]);
  89   parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,NFILE,fnm,
  90                     asize(pa),pa,asize(desc),desc,0,NULL);
  91
  92   top    = read_top(ftp2fn(efTPX,NFILE,fnm));
  93
  94   natoms = read_first_x_v(&status,ftp2fn(efTRN,NFILE,fnm),&t,&x,&v,box);
  95
  96   if (xmax > 0) {
  97     scale  = 5;
  98     nmax   = xmax*scale;
  99   }
 100   else {
 101     scale  = 5;
 102     nmax   = (0.5*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])))*scale;
 103   }
 104   snew(npts,nmax+1);
 105   snew(t_y,nmax+1);
 106   for(i=0; (i<=nmax); i++) {
 107     npts[i] = 0;
 108     t_y[i]  = i/scale;
 109   }
 110   do {
 111     srenew(profile,++nframe);
 112     snew(profile[nframe-1],nmax+1);
 113     srenew(t_x,nframe);
 114     t_x[nframe-1] = t*1000;
 115     hboxx = box[XX][XX]/2;
 116     hboxy = box[YY][YY]/2;
 117     for(i=0; (i<natoms); i++) {
 118       /* determine position dependent on mode */
 119       switch (mode) {
 120       case 0:
 121         xnorm = sqrt(sqr(x[i][XX]-hboxx) + sqr(x[i][YY]-hboxy));
 122         break;
 123       default:
 124         fatal_error(0,"Unknown mode %d",mode);
 125       }
 126       index = xnorm*scale;
 127       if (index <= nmax) {
 128         temp = top->atoms.atom[i].m*iprod(v[i],v[i])/(2*BOLTZ);
 129         if (temp > hi)
 130           hi = temp;
 131         npts[index]++;
 132         profile[nframe-1][index] += temp;
 133       }
 134     }
 135     for(i=0; (i<=nmax); i++) {
 136       if (npts[i] != 0)
 137         profile[nframe-1][i] /= npts[i];
 138       npts[i] = 0;
 139     }
 140   } while (read_next_x_v(status,&t,natoms,x,v,box));
 141   close_trn(status);
 142
 143   fp = ftp2FILE(efXPM,NFILE,fnm,"w");
 144   write_xpm(fp,"Temp. profile","T (a.u.)",
 145             "t (fs)","R (nm)",
 146             nframe,nmax+1,t_x,t_y,profile,0,tmax,
 147             rgblo,rgbhi,&nlevels);
 148
 149   thanx(stderr);
 150
 151   return 0;
 152 }
 153