src/tools/g_rdens.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_rdens_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "sysstuff.h"
  33 #include "string.h"
  34 #include "string2.h"
  35 #include "typedefs.h"
  36 #include "smalloc.h"
  37 #include "macros.h"
  38 #include "princ.h"
  39 #include "vec.h"
  40 #include "xvgr.h"
  41 #include "copyrite.h"
  42 #include "futil.h"
  43 #include "statutil.h"
  44 #include "rdgroup.h"
  45 #include "pbc.h"
  46
  47 real sphere_vol(real r)
  48 {
  49   return (4.0*M_PI/3.0)*(r*r*r);
  50 }
  51
  52 /* rdf_add just counts nr. of atoms in each slice and total mass
  53    in each slice */
  54 void rdf_add(int bin[], real mbin[],real width,real hb2,rvec x[],rvec xcm,
  55              int nx2,atom_id index2[], t_atom atom[])
  56 {
  57   int     j,ind,nout;
  58   int     jx;
  59   rvec    dx;
  60   real    r2;
  61
  62   for(j=0; (j < nx2); j++) {
  63     jx=index2[j];
  64     rvec_sub(x[jx],xcm,dx);
  65     r2=iprod(dx,dx);
  66     if (r2 < hb2) {
  67       ind=sqrt(r2)/width;
  68       mbin[ind]+=atom[jx].m;
  69       bin[ind]++;
  70     }
  71     else
  72       nout++;
  73   }
  74 }
  75
  76 void rdf_calc(char *fn,char *pdens, char *rdens, char *ndens,
  77               real width, int n1,char *grp1,atom_id index1[],
  78               int ngrps,int n2[],char **grpname,atom_id *index[],
  79               t_atom atom[])
  80 {
  81   FILE         *pout,*rout,*nout;
  82   int          **bin;
  83   real         **mbin; /* keeps mass per bin */
  84   char         buf[256];
  85   real         t,hb2,nf_1;
  86   int          i,j,natoms,status,maxbin;
  87   rvec         *x0,xcm;
  88   matrix       box;
  89
  90   if ((natoms=read_first_x(&status,fn,&t,&x0,box))==0)
  91     fatal_error(0,"Could not read coordinates from statusfile\n");
  92
  93   hb2=min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ]))/2;
  94   maxbin=(hb2/width)+2;
  95   snew(bin,ngrps);
  96   snew(mbin,ngrps);
  97   for(i=0; (i<ngrps); i++)
  98   {
  99     snew(bin[i],maxbin);
 100     snew(mbin[i],maxbin);
 101   }
 102
 103   fprintf(stderr,"maxbin = %d, hb2 = %g\n",maxbin,hb2);
 104
 105   j=0;
 106   do {
 107     if ((j % 10) == 0)
 108       fprintf(stderr,"\rframe: %5d",j);
 109
 110     /* Must init pbc every step because of pressure coupling */
 111     init_pbc(box,FALSE);
 112
 113     calc_xcm(x0,n1,index1,atom,xcm,FALSE);
 114
 115     for(i=0; (i<ngrps); i++)
 116       rdf_add(bin[i],mbin[i],width,sqr(hb2),x0,xcm,n2[i],index[i],
 117               atom);
 118     j++;
 119   } while (read_next_x(status,&t,natoms,x0,box));
 120   fprintf(stderr,"\n");
 121   close_trj(status);
 122   sfree(x0);
 123   nf_1=1.0/(real)j;
 124
 125   /* now bin is a mega array with the total nr. of atoms of each type
 126      in each bin. bin*nf_1 is the nr. of atoms of each type in each bin
 127      per frame.
 128      */
 129   for(i=0; (i<ngrps); i++) {
 130     sprintf(buf,"prob:%s-%s.xvg",grp1,grpname[i]);
 131     pout=xvgropen(buf,"Radial Probability Plot","R (nm)","nm\\S-1\\N");
 132     sprintf(buf,"dens:%s-%s.xvg",grp1,grpname[i]);
 133     rout=xvgropen(buf,"Radial Density Plot","R (nm)","g cm\\S-3\\N");
 134     sprintf(buf,"ndens:%s-%s.xvg",grp1,grpname[i]);
 135     nout=xvgropen(buf,"Radial Number Density Plot","R (nm)","Atoms nm\\S-3\\N");
 136
 137     fprintf(pout,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
 138     fprintf(nout,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
 139     fprintf(rout,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
 140
 141     for(j=0; (j<maxbin); j++) {
 142       fprintf(nout,"%10g  %10g\n",width*j,bin[i][j]*nf_1/
 143                 (4*M_PI*width*(j+0.5)*width*(j+0.5)));
 144       fprintf(pout,"%10g  %10g\n",width*j,bin[i][j]*nf_1);
 145       fprintf(rout,"%10g  %10g\n",width*j,mbin[i][j]*nf_1/
 146               (60.22*4*M_PI*width*(j+0.5)*width*(j+0.5)));
 147     }
 148     fclose(nout); fclose(pout); fclose(rout);
 149   }
 150 }
 151
 152 int main(int argc,char *argv[])
 153 {
 154   static char *desc[] = {
 155     "Compute radial densities across the box, in three flavors:"
 156     "probability density, number density, real density"
 157   };
 158   static real      width=0.12;
 159   t_pargs pa[] = {
 160     { "-width", FALSE, etREAL, {&width}, "bin width for radial axis" }
 161   };
 162
 163   char      *grp1,**grpname;                /* groupnames                 */
 164   int       gnx1,*gnx;                      /* sizes of groups            */
 165   atom_id   *ind1,**index;                  /* indices for all groups     */
 166   int       ngrps;
 167   t_topology *top;                          /* topology                   */
 168   t_filenm  fnm[] = {                       /* files for g_order          */
 169     { efTRX, "-f", NULL,  ffREAD },         /* trajectory file            */
 170     { efNDX, NULL, NULL,  ffREAD },         /* index file                 */
 171     { efTPX, NULL, NULL,  ffREAD },         /* topology file              */
 172     { efXVG,"-op","p_rdens", ffWRITE },     /* xvgr output: prob. dens.   */
 173     { efXVG,"-on","n_rdens", ffWRITE },     /* xvgr output: number. dens. */
 174     { efXVG,"-or","r_rdens", ffWRITE },     /* xvgr output: real dens.    */
 175   };
 176
 177 #define NFILE asize(fnm)
 178
 179   CopyRight(stderr,argv[0]);
 180   parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,NFILE,fnm,
 181                     asize(pa),pa,asize(desc),desc,0,NULL);
 182
 183   top = read_top(ftp2fn(efTPX,NFILE,fnm));     /* read topology file */
 184   fprintf(stderr,"Choose first group for Center of Mass computation!\n");
 185
 186   fprintf(stderr,"Select group for Center of Mass calculation:\n");
 187   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx1,&ind1,&grp1);
 188
 189   fprintf(stderr,"How many groups do you want to calc the density plot of?\n");
 190   do {
 191     scanf("%d",&ngrps);
 192   } while (ngrps < 0);
 193
 194   snew(gnx,ngrps);
 195   snew(index,ngrps);
 196   snew(grpname,ngrps);
 197   fprintf(stderr,"Select groups for Density Plot:\n");
 198   rd_index(ftp2fn(efNDX,NFILE,fnm),ngrps,gnx,index,grpname);
 199
 200   rdf_calc(ftp2fn(efTRX,NFILE,fnm),opt2fn("-op",NFILE,fnm),
 201            opt2fn("-or",NFILE,fnm),opt2fn("-on",NFILE,fnm),
 202            width,gnx1,grp1,ind1,ngrps,gnx,grpname,index,top->atoms.atom);
 203
 204   thanx(stdout);
 205
 206   return 0;
 207 }
 208
 209
 210