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