src/contrib/g_anavel.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2006, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "sysstuff.h"
  43 #include "smalloc.h"
  44 #include "macros.h"
  45 #include "statutil.h"
  46 #include "random.h"
  47 #include "names.h"
  48 #include "matio.h"
  49 #include "physics.h"
  50 #include "vec.h"
  51 #include "futil.h"
  52 #include "copyrite.h"
  53 #include "xvgr.h"
  54 #include "string2.h"
  55 #include "index.h"
  56 #include "tpxio.h"
  57
  58 int main(int argc,char *argv[])
  59 {
  60   static char *desc[] = {
  61     "[TT]g_anavel[tt] computes temperature profiles in a sample. The sample",
  62     "can be analysed radial, i.e. the temperature as a function of",
  63     "distance from the center, cylindrical, i.e. as a function of distance",
  64     "from the vector (0,0,1) through the center of the box, or otherwise",
  65     "(will be specified later)"
  66   };
  67   t_filenm fnm[] = {
  68     { efTRN,  "-f",  NULL, ffREAD },
  69     { efTPX,  "-s",  NULL, ffREAD },
  70     { efXPM,  "-o", "xcm", ffWRITE }
  71   };
  72 #define NFILE asize(fnm)
  73
  74   static int  mode = 0,   nlevels = 10;
  75   static real tmax = 300, xmax    = -1;
  76   t_pargs pa[] = {
  77     { "-mode",    FALSE, etINT,  {&mode},    "mode" },
  78     { "-nlevels", FALSE, etINT,  {&nlevels}, "number of levels" },
  79     { "-tmax",    FALSE, etREAL, {&tmax},    "max temperature in output" },
  80     { "-xmax",    FALSE, etREAL, {&xmax},    "max distance from center" }
  81   };
  82
  83   FILE       *fp;
  84   int        *npts,nmax;
  85   int        status;
  86   int        i,j,idum,step,nframe=0,index;
  87   real       temp,rdum,hboxx,hboxy,scale,xnorm=0;
  88   real       **profile=NULL;
  89   real       *t_x=NULL,*t_y,hi=0;
  90   t_topology *top;
  91   int        d,m,n;
  92   matrix     box;
  93   atom_id    *sysindex;
  94   gmx_bool       bHaveV,bReadV;
  95   t_rgb      rgblo = { 0, 0, 1 },rgbhi = { 1, 0, 0 };
  96   int        flags = TRX_READ_X | TRX_READ_V;
  97   t_trxframe fr;
  98
  99
 100   CopyRight(stderr,argv[0]);
 101   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE ,NFILE,fnm,
 102                     asize(pa),pa,asize(desc),desc,0,NULL);
 103
 104   top    = read_top(ftp2fn(efTPX,NFILE,fnm));
 105
 106   read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
 107
 108   if (xmax > 0) {
 109     scale  = 5;
 110     nmax   = xmax*scale;
 111   }
 112   else {
 113     scale  = 5;
 114     nmax   = (0.5*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])))*scale;
 115   }
 116   snew(npts,nmax+1);
 117   snew(t_y,nmax+1);
 118   for(i=0; (i<=nmax); i++) {
 119     npts[i] = 0;
 120     t_y[i]  = i/scale;
 121   }
 122   do {
 123     srenew(profile,++nframe);
 124     snew(profile[nframe-1],nmax+1);
 125     srenew(t_x,nframe);
 126     t_x[nframe-1] = fr.time*1000;
 127     hboxx = box[XX][XX]/2;
 128     hboxy = box[YY][YY]/2;
 129     for(i=0; (i<fr.natoms); i++) {
 130       /* determine position dependent on mode */
 131       switch (mode) {
 132       case 0:
 133         xnorm = sqrt(sqr(fr.x[i][XX]-hboxx) + sqr(fr.x[i][YY]-hboxy));
 134         break;
 135       default:
 136         gmx_fatal(FARGS,"Unknown mode %d",mode);
 137       }
 138       index = xnorm*scale;
 139       if (index <= nmax) {
 140         temp = top->atoms.atom[i].m*iprod(fr.v[i],fr.v[i])/(2*BOLTZ);
 141         if (temp > hi)
 142           hi = temp;
 143         npts[index]++;
 144         profile[nframe-1][index] += temp;
 145       }
 146     }
 147     for(i=0; (i<=nmax); i++) {
 148       if (npts[i] != 0)
 149         profile[nframe-1][i] /= npts[i];
 150       npts[i] = 0;
 151     }
 152   } while (read_next_frame(status,&fr));
 153   close_trx(status);
 154
 155   fp = ftp2FILE(efXPM,NFILE,fnm,"w");
 156   write_xpm(fp,0,"Temp. profile","T (a.u.)",
 157             "t (fs)","R (nm)",
 158             nframe,nmax+1,t_x,t_y,profile,0,tmax,
 159             rgblo,rgbhi,&nlevels);
 160
 161   thanx(stderr);
 162
 163   return 0;
 164 }
 165