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