src/tools/my_rdf.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  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_my_rdf_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "sysstuff.h"
  33 #include "string.h"
  34 #include "typedefs.h"
  35 #include "smalloc.h"
  36 #include "macros.h"
  37 #include "vec.h"
  38 #include "xvgr.h"
  39 #include "pbc.h"
  40 #include "copyrite.h"
  41 #include "futil.h"
  42 #include "statutil.h"
  43 #include "rdgroup.h"
  44
  45 #define RESOLUTION 0.01   /* 0.01 nm resolution */
  46 #define NR_HIST 200      /* compile histogram of 200 bins (from 0 to 2.0 nm)*/
  47
  48 static real hist[NR_HIST];
  49 static real dens_tot;
  50
  51 /****************************************************************************/
  52 /* This program calculates density corrected and uncorrected radial         */
  53 /* distribution functions in inhomogeneous systems like membranes.          */
  54 /* Peter Tieleman, January 1996                                             */
  55 /****************************************************************************/
  56
  57 real sphere_vol(real r)
  58 {
  59   return (4.0*M_PI/3.0)*(r*r*r);
  60 }
  61
  62 void calc_rdf(char *fn, atom_id **index, int gnx[], real l)
  63 {
  64   rvec *x0;                /* coordinates without pbc */
  65   matrix box;              /* box (3x3) */
  66   int natoms,              /* nr. atoms in trj */
  67     status,
  68     i,j,m,                 /* loop indices */
  69     teller = 0,
  70     nr_frames = 0,         /* number of frames */
  71     nwater, nlipid;        /* number of waters, number of lipids */
  72   real t,
  73     rho_local;
  74   int bin;                 /* bin to put atom in */
  75   atom_id *water, *lipid;  /* the index numbers for lipid and water */
  76   rvec tmp;
  77   real dist;
  78   real count = 0; real dens = 0;
  79
  80   for (i = 0; i < NR_HIST; i++)
  81     hist[i] = 0;
  82
  83   if ((natoms = read_first_x(&status,fn,&t,&x0,box)) == 0)
  84     fatal_error(0,"Could not read coordinates from statusfile\n");
  85
  86   fprintf(stderr,"Cut-off for counting is %f\n",l);
  87   nwater = gnx[1]; nlipid = gnx[0];
  88
  89   rho_local = nwater/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
  90   fprintf(stderr,"Average water density: %8.3g\n"
  91           "Number of lipids: %d. Number of waters: %d\n",
  92           rho_local, nlipid, nwater);
  93
  94   lipid = index[0];
  95   water = index[1];
  96
  97   dens_tot = 0;
  98
  99   /*********** Start processing trajectory ***********/
 100   do
 101   {
 102     if ((teller++ % 10) == 0)
 103       fprintf(stderr,"\rFrame: %d",teller-1);
 104
 105     for (i = 0; i < nlipid; i++)
 106     {
 107       for (j = 0; j < nwater; j++)  /* nice double loop */
 108       {
 109         rvec_sub(x0[lipid[i]], x0[water[j]], tmp);
 110         for (m = 0; m < DIM; m++)
 111         {
 112           if (tmp[m] > 0.5*box[m][m])
 113             tmp[m] -= box[m][m];
 114           else if (tmp[m] < -0.5*box[m][m])
 115             tmp[m] += box[m][m];
 116         }
 117         dist = norm(tmp);
 118         /*      fprintf(stderr,"dist: %f",dist); */
 119         if (dist < RESOLUTION*NR_HIST)
 120         {
 121           bin = trunc(dist/RESOLUTION);
 122           if (bin < 0 || bin > NR_HIST)
 123           {
 124             fprintf(stderr,"Warning: bin out of range: %d\n",bin);
 125             bin = NR_HIST;
 126           }
 127           hist[bin] += 1.0;
 128         }
 129         if (dist < 2.0)
 130           dens += 1.0;
 131
 132         if (dist < l)
 133         {
 134           count += 1.0;
 135         }
 136       }
 137       dens_tot += dens;
 138       dens = 0;
 139     }
 140   } while (read_next_x(status,&t,natoms,x0,box));
 141
 142   /*********** done with status file **********/
 143   close_trj(status);
 144   dens_tot /= nlipid * teller * sphere_vol(2);
 145
 146   fprintf(stderr,"Local density of water around lipid: %f\n",dens_tot);
 147
 148   for(i=1;i<NR_HIST;i++)
 149     hist[i] /= nlipid*teller*dens_tot;
 150
 151   count /= nlipid*teller;
 152   fprintf(stderr,"Counted %g OW's\n", count);
 153   sfree(x0);  /* free memory used by coordinate array */
 154 }
 155
 156 real *integrate(real data[])
 157 {
 158   int i,j;
 159   real *result,
 160     sum;
 161
 162   snew(result, NR_HIST);
 163
 164   for (i = 1; i < NR_HIST ; i++)
 165   {
 166     sum = 0;
 167     for (j = 1; j < i; j++)
 168       sum += data[j] + 0.5 * (data[j+1] - data[j]);
 169     result[i] = sum * dens_tot;
 170   }
 171   return result;
 172 }
 173
 174 void plot_rdf(char *afile, char *grpname[])
 175 {
 176   FILE       *uncor;     /* xvgr file corrected rdf   */
 177   char       buf[256]; /* for xvgr title */
 178   int        i, n;
 179   real       *integ;
 180
 181   sprintf(buf,"Uncorrected rdf");
 182   uncor = xvgropen(afile, buf, "r", "g(r)");
 183
 184   integ = integrate(hist);
 185
 186   for (i = 0; i < NR_HIST; i++)
 187   {
 188     hist[i]/= 4*M_PI*RESOLUTION*RESOLUTION*RESOLUTION*i*i;
 189   /* r2 = RESOLUTION*RESOLUTION*i*i. The third comes from the change of the x-ax
 190      from i to i*RESOLUTION. sucks */
 191     fprintf(uncor,"%12g   %12g   %12g\n", i*RESOLUTION, hist[i], integ[i]);
 192   }
 193
 194   fclose(uncor);
 195 }
 196
 197 void main(int argc,char *argv[])
 198 {
 199   static char *desc[] = {
 200     "Compute rdf's"
 201   };
 202   static char *opts[] = {
 203     "-l distance"
 204   };
 205   static char *odesc[] = {
 206     "distance taken as minimum for counting"
 207   };
 208
 209   t_manual man = {asize(desc),desc,asize(opts),opts,NULL,NULL};
 210   real      l = 0;
 211   char      **grpname;                      /* groupnames                 */
 212   int       i,ngrps = 0,                    /* nr. of groups              */
 213     *ngx;                                   /* sizes of groups            */
 214   atom_id   **index;                        /* indices for all groups     */
 215   t_filenm  fnm[] = {                       /* files for g_order          */
 216     { efTRX, "-f", NULL,  ffREAD },         /* trajectory file            */
 217     { efNDX, NULL, NULL,  ffREAD },         /* index file                 */
 218     { efXVG, "-o", "graph",ffWRITE },       /* xvgr output file           */
 219   };
 220
 221 #define NFILE asize(fnm)
 222
 223   CopyRight(stderr,argv[0]);
 224   parse_common_args(&argc, argv, PCA_CAN_TIME, NFILE, fnm, TRUE, &man);
 225
 226   for (i = 1; i < argc; i++) {
 227     if (strcmp(argv[i],"-l") == 0)  {
 228       if (i < argc - 1) {
 229         sscanf(argv[i+1],"%f",&l);
 230         i++;
 231       }
 232     }
 233   }
 234
 235   snew(grpname,2);
 236   snew(index,2);
 237   snew(ngx,2);
 238
 239   rd_index(ftp2fn(efNDX,NFILE,fnm),2,ngx,index,grpname);
 240
 241   fprintf(stderr,"Assuming group %s is the solvent.\n", grpname[1]);
 242
 243   calc_rdf(ftp2fn(efTRX,NFILE,fnm),index, ngx, l);
 244
 245   plot_rdf(opt2fn("-o",NFILE,fnm), grpname);
 246 }