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