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