3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Groningen Machine for Chemical Simulation
39 #include "gromacs/utility/smalloc.h"
41 #include "gromacs/commandline/pargs.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"
49 #include "gromacs/fileio/tpxio.h"
51 int main(int argc,char *argv[])
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)"
61 { efTRN, "-f", NULL, ffREAD },
62 { efTPX, "-s", NULL, ffREAD },
63 { efXPM, "-o", "xcm", ffWRITE }
65 #define NFILE asize(fnm)
67 static int mode = 0, nlevels = 10;
68 static real tmax = 300, xmax = -1;
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" }
79 int i,j,idum,step,nframe=0,index;
80 real temp,rdum,hboxx,hboxy,scale,xnorm=0;
82 real *t_x=NULL,*t_y,hi=0;
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;
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);
97 top = read_top(ftp2fn(efTPX,NFILE,fnm));
99 read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
107 nmax = (0.5*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])))*scale;
111 for(i=0; (i<=nmax); i++) {
116 srenew(profile,++nframe);
117 snew(profile[nframe-1],nmax+1);
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 */
126 xnorm = sqrt(sqr(fr.x[i][XX]-hboxx) + sqr(fr.x[i][YY]-hboxy));
129 gmx_fatal(FARGS,"Unknown mode %d",mode);
133 temp = top->atoms.atom[i].m*iprod(fr.v[i],fr.v[i])/(2*BOLTZ);
137 profile[nframe-1][index] += temp;
140 for(i=0; (i<=nmax); i++) {
142 profile[nframe-1][i] /= npts[i];
145 } while (read_next_frame(status,&fr));
148 fp = ftp2FILE(efXPM,NFILE,fnm,"w");
149 write_xpm(fp,0,"Temp. profile","T (a.u.)",
151 nframe,nmax+1,t_x,t_y,profile,0,tmax,
152 rgblo,rgbhi,&nlevels);