src/tools/g_clustsize.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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33 #include <math.h>
  34 #include <ctype.h>
  35 #include "string2.h"
  36 #include "sysstuff.h"
  37 #include "typedefs.h"
  38 #include "macros.h"
  39 #include "vec.h"
  40 #include "pbc.h"
  41 #include "rmpbc.h"
  42 #include "statutil.h"
  43 #include "xvgr.h"
  44 #include "copyrite.h"
  45 #include "futil.h"
  46 #include "statutil.h"
  47 #include "tpxio.h"
  48 #include "rdgroup.h"
  49 #include "smalloc.h"
  50 #include "fftgrid.h"
  51 #include "calcgrid.h"
  52 #include "nrnb.h"
  53 #include "shift_util.h"
  54 #include "pme.h"
  55 #include "gstat.h"
  56 #include "matio.h"
  57
  58 static void clust_size(char *ndx,char *trx,char *xpm,char *ncl,char *acl,
  59                        char *mcl,real cut,int nskip,int nlevels,
  60                        t_rgb rmid,t_rgb rhi)
  61 {
  62   FILE    *fp,*gp,*hp;
  63   atom_id *index=NULL;
  64   int     nindex,natoms,status;
  65   rvec    *x=NULL,dx;
  66   matrix  box;
  67   char    *gname;
  68   real    t,dx2,cut2,**cs_dist=NULL,*t_x=NULL,*t_y,mid,cav;
  69   int     i,j,k,ai,aj,ak,ci,cj,nframe,nclust,n_x,n_y,max_size=0;
  70   int     *clust_index,*clust_size,max_clust_size,nav;
  71   t_rgb   rlo = { 1.0, 1.0, 1.0 };
  72
  73   fp = xvgropen(ncl,"Number of clusters","Time (ps)","N");
  74   gp = xvgropen(acl,"Average cluster size","Time (ps)","#molecules");
  75   hp = xvgropen(mcl,"Max cluster size","Time (ps)","#molecules");
  76   rd_index(ndx,1,&nindex,&index,&gname);
  77   natoms = read_first_x(&status,trx,&t,&x,box);
  78   snew(clust_index,nindex);
  79   snew(clust_size,nindex);
  80   cut2   = cut*cut;
  81   nframe = 0;
  82   n_x    = 0;
  83   snew(t_y,nindex);
  84   for(i=0; (i<nindex); i++)
  85     t_y[i] = i+1;
  86   do {
  87     if ((nskip == 0) || ((nskip > 0) && ((nframe % nskip) == 0))) {
  88       init_pbc(box);
  89       max_clust_size = 1;
  90       for(i=0; (i<nindex); i++) {
  91         clust_index[i] = i;
  92         clust_size[i]  = 1;
  93       }
  94
  95       for(i=0; (i<nindex); i++) {
  96         ai = index[i];
  97         ci = clust_index[i];
  98         for(j=i+1; (j<nindex); j++) {
  99           cj = clust_index[j];
 100           if (ci != cj) {
 101             aj = index[j];
 102             pbc_dx(x[ai],x[aj],dx);
 103             dx2 = iprod(dx,dx);
 104             if (dx2 < cut2) {
 105               /* Merge clusters */
 106               for(k=j; (k<nindex); k++) {
 107                 if (clust_index[k] == cj) {
 108                   clust_size[cj]--;
 109                   clust_index[k] = ci;
 110                   clust_size[i]++;
 111                 }
 112               }
 113             }
 114           }
 115         }
 116       }
 117       n_x++;
 118       srenew(t_x,n_x);
 119       t_x[n_x-1] = t;
 120       srenew(cs_dist,n_x);
 121       snew(cs_dist[n_x-1],nindex);
 122       nclust = 0;
 123       cav    = 0;
 124       nav    = 0;
 125       for(i=0; (i<nindex); i++) {
 126         ci = clust_size[i];
 127         if (ci > max_clust_size) max_clust_size = ci;
 128         if (ci > 0) {
 129           nclust++;
 130           cs_dist[n_x-1][ci-1] += 100.0*ci/(real)nindex;
 131           max_size = max(max_size,ci);
 132           if (ci > 1) {
 133             cav += ci;
 134             nav++;
 135           }
 136         }
 137       }
 138       fprintf(fp,"%10.3e  %10d\n",t,nclust);
 139       if (nav > 0)
 140         fprintf(gp,"%10.3e  %10.3f\n",t,cav/nav);
 141       fprintf(hp, "%10.3e  %10d\n", t, max_clust_size);
 142     }
 143     nframe++;
 144   } while (read_next_x(status,&t,natoms,x,box));
 145   close_trx(status);
 146   fclose(fp);
 147   fclose(gp);
 148   fclose(hp);
 149
 150   /* Look for the smallest entry that is not zero */
 151   mid = 100.0;
 152   for(i=0; (i<n_x); i++)
 153     for(j=0; (j<max_size); j++)
 154       if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < mid))
 155         mid = cs_dist[i][j];
 156
 157   fp = ffopen(xpm,"w");
 158   write_xpm3(fp,"Cluster size distribution","Fraction","Time (ps)","Size",
 159              n_x,max_size-1,t_x,t_y,cs_dist,0,mid,100.0,
 160              rlo,rmid,rhi,&nlevels);
 161   fclose(fp);
 162
 163   sfree(clust_index);
 164   sfree(clust_size);
 165   sfree(index);
 166 }
 167
 168 int main(int argc,char *argv[])
 169 {
 170   static char *desc[] = {
 171     "This program computes the size distributions of molecular/atomic clusters in",
 172     "the gas phase. The output is given in the form of a XPM file.",
 173     "The total number of clusters is written to a XVG file."
 174   };
 175   static real cutoff   = 0.35;
 176   static int  nskip    = 0;
 177   static int  nlevels  = 20;
 178   static rvec rlo      = { 1.0, 1.0, 0.0 };
 179   static rvec rhi      = { 0.0, 0.0, 1.0 };
 180   t_pargs pa[] = {
 181     { "-cut",      FALSE, etREAL, {&cutoff},
 182       "Largest distance (nm) to be considered in a cluster" },
 183     { "-nskip",    FALSE, etINT,  {&nskip},
 184       "Number of frames to skip between writing" },
 185     { "-nlevels",  FALSE, etINT,  {&nlevels},
 186       "Number of levels of grey in xpm output" },
 187     { "-rgblo",    FALSE, etRVEC, {rlo},
 188       "RGB values for the color of the lowest occupied cluster size" },
 189     { "-rgbhi",    FALSE, etRVEC, {rhi},
 190       "RGB values for the color of the highest occupied cluster size" }
 191   };
 192 #define NPA asize(pa)
 193   char       *fnTPS,*fnNDX;
 194   bool       bSQ,bRDF;
 195   t_rgb      rgblo,rgbhi;
 196
 197   t_filenm   fnm[] = {
 198     { efTRX, "-f",  NULL,    ffREAD  },
 199     { efNDX, NULL,  NULL,    ffOPTRD },
 200     { efXPM, "-o", "csize",  ffWRITE },
 201     { efXVG, "-nc","nclust", ffWRITE },
 202     { efXVG, "-mc","maxclust", ffWRITE },
 203     { efXVG, "-ac","avclust", ffWRITE }
 204   };
 205 #define NFILE asize(fnm)
 206
 207   CopyRight(stderr,argv[0]);
 208   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 209                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
 210
 211   fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
 212   rgblo.r = rlo[XX],rgblo.g = rlo[YY],rgblo.b = rlo[ZZ];
 213   rgbhi.r = rhi[XX],rgbhi.g = rhi[YY],rgbhi.b = rhi[ZZ];
 214
 215   clust_size(fnNDX,ftp2fn(efTRX,NFILE,fnm),ftp2fn(efXPM,NFILE,fnm),
 216              opt2fn("-nc",NFILE,fnm),opt2fn("-ac",NFILE,fnm),opt2fn("-mc",NFILE,fnm),
 217              cutoff,nskip,nlevels,rgblo,rgbhi);
 218
 219   thanx(stderr);
 220
 221   return 0;
 222 }