src/tools/g_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  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_rdf_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <ctype.h>
  33 #include "string2.h"
  34 #include "sysstuff.h"
  35 #include "typedefs.h"
  36 #include "macros.h"
  37 #include "vec.h"
  38 #include "pbc.h"
  39 #include "rmpbc.h"
  40 #include "xvgr.h"
  41 #include "copyrite.h"
  42 #include "futil.h"
  43 #include "statutil.h"
  44 #include "tpxio.h"
  45 #include "rdgroup.h"
  46 #include "smalloc.h"
  47
  48 int main(int argc,char *argv[])
  49 {
  50   static char *desc[] = {
  51     "g_rdf calculates radial distribution functions in different ways.",
  52     "The normal method is around a (set of) particle(s), the other method",
  53     "is around the center of mass of a set of particles.[PAR]",
  54     "If a run input file is supplied ([TT]-s[tt]), exclusions defined",
  55     "in that file are taken into account when calculating the rdf.",
  56     "The option [TT]-cut[tt] is meant as an alternative way to avoid",
  57     "intramolecular peaks in the rdf plot.",
  58     "It is however better to supply a run input file with a higher number of",
  59     "exclusions. For eg. benzene a topology with nrexcl set to 5",
  60     "would eliminate all intramolecular contributions to the rdf.",
  61     "Note that all atoms in the selected groups are used, also the ones",
  62     "that don't have Lennard-Jones interactions."
  63   };
  64   static bool bCM=FALSE;
  65   static real cutoff=0, binwidth=0.005;
  66   t_pargs pa[] = {
  67     { "-bin",      FALSE, etREAL, {&binwidth},
  68       "Binwidth (nm)" },
  69     { "-com",      FALSE, etBOOL, {&bCM},
  70       "RDF with respect to the center of mass of first group" },
  71     { "-cut",      FALSE, etREAL, {&cutoff},
  72       "Shortest distance (nm) to be considered"},
  73   };
  74 #define NPA asize(pa)
  75   FILE       *fp;
  76   char       *fnTPS,*fnNDX;
  77   int        status;
  78   char       **grpname;
  79   char       outf1[STRLEN],outf2[STRLEN];
  80   char       title[STRLEN];
  81   int        natoms,i,j,k,nbin,j0,j1,n;
  82   int        *count;
  83   unsigned long int sum;
  84   real       t,boxmin,hbox,hbox2,cut2,r,r2,invbinw,normfac;
  85   real       segvol,spherevol,prev_spherevol;
  86   rvec       *x,xcom,dx;
  87   real       *inv_segvol;
  88   bool       *bExcl,bTop;
  89   matrix     box;
  90   int        isize[3],*npairs;
  91   atom_id    *index[3],ix,jx,**pairs;
  92   t_topology top;
  93   t_block    *excl;
  94   t_filenm   fnm[] = {
  95     { efTRX, "-f", NULL,   ffREAD },
  96     { efTPS, NULL, NULL,   ffOPTRD },
  97     { efNDX, NULL, NULL,   ffOPTRD },
  98     { efXVG, NULL, "rdf",  ffWRITE },
  99   };
 100 #define NFILE asize(fnm)
 101
 102   CopyRight(stderr,argv[0]);
 103   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 104                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
 105
 106   fnTPS = ftp2fn_null(efTPS,NFILE,fnm);
 107   fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
 108   if (!fnTPS && !fnNDX)
 109     fatal_error(0,"Neither index file nor topology file specified\n"
 110                 "             Nothing to do!");
 111
 112   excl=NULL;
 113   if (fnTPS) {
 114     bTop=read_tps_conf(fnTPS,title,&top,&x,NULL,box,TRUE);
 115     mk_single_top(&top);
 116     if (bTop && !bCM)
 117       /* get exclusions from topology */
 118       excl=&(top.atoms.excl);
 119   }
 120
 121   snew(grpname,2);
 122   if (fnTPS)
 123     get_index(&top.atoms,fnNDX,2,isize,index,grpname);
 124   else
 125     rd_index(fnNDX,2,isize,index,grpname);
 126
 127   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 128   if ( !natoms )
 129     fatal_error(0,"Could not read coordinates from statusfile\n");
 130   if (fnTPS)
 131     /* check with topology */
 132     if ( natoms > top.atoms.nr )
 133       fatal_error(0,"Trajectory (%d atoms) does not match topology (%d atoms)",
 134                   natoms,top.atoms.nr);
 135   /* check with index groups */
 136   for (i=0; i<2; i++)
 137     for (j=0; j<isize[i]; j++)
 138       if ( index[i][j] >= natoms )
 139         fatal_error(0,"Atom index (%d) in index group %s (%d atoms) larger "
 140                     "than number of atoms in trajectory (%d atoms)",
 141                     index[i][j],grpname[i],isize[i],natoms);
 142
 143   if (bCM) {
 144     /* move index[0] to index[2] and make 'dummy' index[0] */
 145     isize[2]=isize[0];
 146     snew(index[2],isize[2]);
 147     for(i=0; i<isize[0]; i++)
 148       index[2][i]=index[0][i];
 149     isize[0]=1;
 150     index[0][0]=natoms;
 151     srenew(index[0],isize[0]);
 152     /* make space for center of mass */
 153     srenew(x,natoms+1);
 154   }
 155
 156   /* initialize some handy things */
 157   boxmin = min( norm(box[XX]), min( norm(box[YY]), norm(box[ZZ]) ) );
 158   hbox   = boxmin / 2.0;
 159   nbin   = (int)(hbox / binwidth) + 1;
 160   invbinw = 1.0 / binwidth;
 161   hbox2  = sqr(hbox);
 162   cut2   = sqr(cutoff);
 163
 164   /* this is THE array */
 165   snew(count,nbin+1);
 166
 167   /* make pairlist array for groups and exclusions */
 168   snew(pairs,isize[0]);
 169   snew(npairs,isize[0]);
 170   for(i=0; i<isize[0]; i++)
 171     snew(pairs[i], isize[1]);
 172   snew(bExcl,natoms);
 173   for( i = 0; i < isize[0]; i++) {
 174     ix = index[0][i];
 175     for( j = 0; j < natoms; j++)
 176       bExcl[j] = FALSE;
 177     /* exclusions? */
 178     if (excl)
 179       for( j = excl->index[ix]; j < excl->index[ix+1]; j++)
 180         bExcl[excl->a[j]]=TRUE;
 181     k = 0;
 182     for( j = 0; j < isize[1]; j++) {
 183       jx = index[1][j];
 184       if (!bExcl[jx])
 185         pairs[i][k++]=jx;
 186     }
 187     npairs[i]=k;
 188     srenew(pairs[i],npairs[i]);
 189   }
 190   sfree(bExcl);
 191
 192   do {
 193     /* Must init pbc every step because of pressure coupling */
 194     init_pbc(box,FALSE);
 195     rm_pbc(&top.idef,natoms,box,x,x);
 196
 197     if (bCM) {
 198       /* calculate centre of mass */
 199       clear_rvec(xcom);
 200       for(i=0; (i < isize[2]); i++) {
 201         ix = index[2][i];
 202         rvec_inc(xcom,x[ix]);
 203       }
 204       /* store it in the first 'group' */
 205       for(j=0; (j<DIM); j++)
 206         x[index[0][0]][j] = xcom[j] / isize[2];
 207     }
 208
 209     for(i=0; i < isize[0]; i++) {
 210       ix = index[0][i];
 211       for(j=0; j < npairs[i]; j++) {
 212         jx=pairs[i][j];
 213         pbc_dx(x[ix],x[jx],dx);
 214         r2=iprod(dx,dx);
 215         if ( r2 <= hbox2 && r2 > cut2 )
 216           count[(int)(sqrt(r2)*invbinw)]++;
 217       }
 218     }
 219   } while (read_next_x(status,&t,natoms,x,box));
 220   fprintf(stderr,"\n");
 221
 222   close_trj(status);
 223
 224   sfree(x);
 225
 226   /* Calculate volume of sphere segments */
 227   snew(inv_segvol,nbin);
 228   prev_spherevol=0;
 229   for(i=0; (i<nbin); i++) {
 230     r = (i+1)*binwidth;
 231     spherevol=(4.0/3.0)*M_PI*r*r*r;
 232     segvol=spherevol-prev_spherevol;
 233     inv_segvol[i]=1.0/segvol;
 234     prev_spherevol=spherevol;
 235   }
 236
 237   /* Calculate normalization factor and totals */
 238   sum = 0;
 239   for(i=0; i<nbin; i++)
 240     sum += count[i];
 241   r = nbin*binwidth;
 242   normfac = (4.0/3.0)*M_PI * r*r*r / sum;
 243
 244   fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),"Radial Distribution","r","");
 245   fprintf(fp,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
 246   sum = 0;
 247   for(i = 0; i < nbin; i++) {
 248     sum += count[i];
 249     fprintf(fp, "%10g %10g %10g\n", (i+0.5)*binwidth,
 250             count[i]*inv_segvol[i]*normfac, (real)sum*normfac);
 251   }
 252   ffclose(fp);
 253
 254   xvgr_file(ftp2fn(efXVG,NFILE,fnm),NULL);
 255
 256   thanx(stdout);
 257
 258   return 0;
 259 }