src/tools/gmx_clustsize.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.2
  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-2007, 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 <math.h>
  40 #include <ctype.h>
  41
  42 #include "string2.h"
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "macros.h"
  46 #include "vec.h"
  47 #include "pbc.h"
  48 #include "rmpbc.h"
  49 #include "statutil.h"
  50 #include "xvgr.h"
  51 #include "copyrite.h"
  52 #include "futil.h"
  53 #include "statutil.h"
  54 #include "tpxio.h"
  55 #include "index.h"
  56 #include "smalloc.h"
  57 #include "calcgrid.h"
  58 #include "nrnb.h"
  59 #include "physics.h"
  60 #include "coulomb.h"
  61 #include "pme.h"
  62 #include "gstat.h"
  63 #include "matio.h"
  64 #include "mtop_util.h"
  65 #include "gmx_ana.h"
  66
  67
  68 static void clust_size(const char *ndx,const char *trx,const char *xpm,
  69                        const char *xpmw,const char *ncl,const char *acl,
  70                        const char *mcl,const char *histo,const char *tempf,
  71                        const char *mcn,bool bMol,bool bPBC,const char *tpr,
  72                        real cut,int nskip,int nlevels,
  73                        t_rgb rmid,t_rgb rhi,int ndf,
  74                        const output_env_t oenv)
  75 {
  76   FILE    *fp,*gp,*hp,*tp;
  77   atom_id *index=NULL;
  78   int     nindex,natoms,status;
  79   rvec    *x=NULL,*v=NULL,dx;
  80   t_pbc   pbc;
  81   char    *gname;
  82   char    timebuf[32];
  83   bool    bSame,bTPRwarn=TRUE;
  84   /* Topology stuff */
  85   t_trxframe  fr;
  86   t_tpxheader tpxh;
  87   gmx_mtop_t *mtop=NULL;
  88   int     ePBC=-1;
  89   t_block *mols=NULL;
  90   t_atom  *atom;
  91   int     version,generation,ii,jj,nsame;
  92   real    temp,tfac;
  93   /* Cluster size distribution (matrix) */
  94   real    **cs_dist=NULL;
  95   real    tf,dx2,cut2,*t_x=NULL,*t_y,cmid,cmax,cav,ekin;
  96   int     i,j,k,ai,aj,ak,ci,cj,nframe,nclust,n_x,n_y,max_size=0;
  97   int     *clust_index,*clust_size,max_clust_size,max_clust_ind,nav,nhisto;
  98   t_rgb   rlo = { 1.0, 1.0, 1.0 };
  99
 100   clear_trxframe(&fr,TRUE);
 101   sprintf(timebuf,"Time (%s)",output_env_get_time_unit(oenv));
 102   tf     = output_env_get_time_factor(oenv);
 103   fp     = xvgropen(ncl,"Number of clusters",timebuf,"N",oenv);
 104   gp     = xvgropen(acl,"Average cluster size",timebuf,"#molecules",oenv);
 105   hp     = xvgropen(mcl,"Max cluster size",timebuf,"#molecules",oenv);
 106   tp     = xvgropen(tempf,"Temperature of largest cluster",timebuf,"T (K)",
 107                     oenv);
 108
 109   if (!read_first_frame(oenv,&status,trx,&fr,TRX_NEED_X | TRX_READ_V))
 110     gmx_file(trx);
 111
 112   natoms = fr.natoms;
 113   x      = fr.x;
 114
 115   if (tpr) {
 116     snew(mtop,1);
 117     read_tpxheader(tpr,&tpxh,TRUE,&version,&generation);
 118     if (tpxh.natoms != natoms)
 119       gmx_fatal(FARGS,"tpr (%d atoms) and xtc (%d atoms) do not match!",
 120                 tpxh.natoms,natoms);
 121     ePBC = read_tpx(tpr,NULL,NULL,&natoms,NULL,NULL,NULL,mtop);
 122   }
 123   if (ndf <= -1)
 124     tfac = 1;
 125   else
 126     tfac = ndf/(3.0*natoms);
 127
 128   if (bMol) {
 129     if (ndx)
 130       printf("Using molecules rather than atoms. Not reading index file %s\n",
 131              ndx);
 132     mols = &(mtop->mols);
 133
 134     /* Make dummy index */
 135     nindex = mols->nr;
 136     snew(index,nindex);
 137     for(i=0; (i<nindex); i++)
 138       index[i] = i;
 139     gname = strdup("mols");
 140   }
 141   else
 142     rd_index(ndx,1,&nindex,&index,&gname);
 143
 144   snew(clust_index,nindex);
 145   snew(clust_size,nindex);
 146   cut2   = cut*cut;
 147   nframe = 0;
 148   n_x    = 0;
 149   snew(t_y,nindex);
 150   for(i=0; (i<nindex); i++)
 151     t_y[i] = i+1;
 152   max_clust_size = 1;
 153   max_clust_ind  = -1;
 154   do {
 155     if ((nskip == 0) || ((nskip > 0) && ((nframe % nskip) == 0))) {
 156       if (bPBC)
 157         set_pbc(&pbc,ePBC,fr.box);
 158       max_clust_size = 1;
 159       max_clust_ind  = -1;
 160
 161       /* Put all atoms/molecules in their own cluster, with size 1 */
 162       for(i=0; (i<nindex); i++) {
 163         /* Cluster index is indexed with atom index number */
 164         clust_index[i] = i;
 165         /* Cluster size is indexed with cluster number */
 166         clust_size[i]  = 1;
 167       }
 168
 169       /* Loop over atoms */
 170       for(i=0; (i<nindex); i++) {
 171         ai = index[i];
 172         ci = clust_index[i];
 173
 174         /* Loop over atoms (only half a matrix) */
 175         for(j=i+1; (j<nindex); j++) {
 176           cj = clust_index[j];
 177
 178           /* If they are not in the same cluster already */
 179           if (ci != cj) {
 180             aj = index[j];
 181
 182             /* Compute distance */
 183             if (bMol) {
 184               bSame = FALSE;
 185               for(ii=mols->index[ai]; !bSame && (ii<mols->index[ai+1]); ii++) {
 186                 for(jj=mols->index[aj]; !bSame && (jj<mols->index[aj+1]); jj++) {
 187                   if (bPBC)
 188                     pbc_dx(&pbc,x[ii],x[jj],dx);
 189                   else
 190                     rvec_sub(x[ii],x[jj],dx);
 191                   dx2   = iprod(dx,dx);
 192                   bSame = (dx2 < cut2);
 193                 }
 194               }
 195             }
 196             else {
 197               if (bPBC)
 198                 pbc_dx(&pbc,x[ai],x[aj],dx);
 199               else
 200                 rvec_sub(x[ai],x[aj],dx);
 201               dx2 = iprod(dx,dx);
 202               bSame = (dx2 < cut2);
 203             }
 204             /* If distance less than cut-off */
 205             if (bSame) {
 206               /* Merge clusters: check for all atoms whether they are in
 207                * cluster cj and if so, put them in ci
 208                */
 209               for(k=0; (k<nindex); k++) {
 210                 if ((clust_index[k] == cj)) {
 211                   if (clust_size[cj] <= 0)
 212                     gmx_fatal(FARGS,"negative cluster size %d for element %d",
 213                                 clust_size[cj],cj);
 214                   clust_size[cj]--;
 215                   clust_index[k] = ci;
 216                   clust_size[ci]++;
 217                 }
 218               }
 219             }
 220           }
 221         }
 222       }
 223       n_x++;
 224       srenew(t_x,n_x);
 225       t_x[n_x-1] = fr.time*tf;
 226       srenew(cs_dist,n_x);
 227       snew(cs_dist[n_x-1],nindex);
 228       nclust = 0;
 229       cav    = 0;
 230       nav    = 0;
 231       for(i=0; (i<nindex); i++) {
 232         ci = clust_size[i];
 233         if (ci > max_clust_size) {
 234           max_clust_size = ci;
 235           max_clust_ind  = i;
 236         }
 237         if (ci > 0) {
 238           nclust++;
 239           cs_dist[n_x-1][ci-1] += 1.0;
 240           max_size = max(max_size,ci);
 241           if (ci > 1) {
 242             cav += ci;
 243             nav++;
 244           }
 245         }
 246       }
 247       fprintf(fp,"%14.6e  %10d\n",fr.time,nclust);
 248       if (nav > 0)
 249         fprintf(gp,"%14.6e  %10.3f\n",fr.time,cav/nav);
 250       fprintf(hp, "%14.6e  %10d\n",fr.time,max_clust_size);
 251     }
 252     /* Analyse velocities, if present */
 253     if (fr.bV) {
 254       if (!tpr) {
 255         if (bTPRwarn) {
 256           printf("You need a tpr file to analyse temperatures\n");
 257           bTPRwarn = FALSE;
 258         }
 259       }
 260       else {
 261         v = fr.v;
 262         /* Loop over clusters and for each cluster compute 1/2 m v^2 */
 263         if (max_clust_ind >= 0) {
 264           ekin = 0;
 265           for(i=0; (i<nindex); i++)
 266             if (clust_index[i] == max_clust_ind) {
 267               ai    = index[i];
 268               gmx_mtop_atomnr_to_atom(mtop,ai,&atom);
 269               ekin += 0.5*atom->m*iprod(v[ai],v[ai]);
 270             }
 271           temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
 272           fprintf(tp,"%10.3f  %10.3f\n",fr.time,temp);
 273         }
 274       }
 275     }
 276     nframe++;
 277   } while (read_next_frame(oenv,status,&fr));
 278   close_trx(status);
 279   ffclose(fp);
 280   ffclose(gp);
 281   ffclose(hp);
 282   ffclose(tp);
 283   if (max_clust_ind >= 0) {
 284     fp = ffopen(mcn,"w");
 285     fprintf(fp,"[ max_clust ]\n");
 286     for(i=0; (i<nindex); i++)
 287       if (clust_index[i] == max_clust_ind) {
 288         if (bMol) {
 289           for(j=mols->index[i]; (j<mols->index[i+1]); j++)
 290             fprintf(fp,"%d\n",j+1);
 291         }
 292         else {
 293           fprintf(fp,"%d\n",index[i]+1);
 294         }
 295       }
 296     ffclose(fp);
 297   }
 298
 299   /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
 300   fp = xvgropen(histo,"Cluster size distribution","Cluster size","()",oenv);
 301   nhisto = 0;
 302   fprintf(fp,"%5d  %8.3f\n",0,0.0);
 303   for(j=0; (j<max_size); j++) {
 304     real nelem = 0;
 305     for(i=0; (i<n_x); i++)
 306       nelem += cs_dist[i][j];
 307     fprintf(fp,"%5d  %8.3f\n",j+1,nelem/n_x);
 308     nhisto += (int)((j+1)*nelem/n_x);
 309   }
 310   fprintf(fp,"%5d  %8.3f\n",j+1,0.0);
 311   ffclose(fp);
 312
 313   fprintf(stderr,"Total number of atoms in clusters =  %d\n",nhisto);
 314
 315   /* Look for the smallest entry that is not zero
 316    * This will make that zero is white, and not zero is coloured.
 317    */
 318   cmid = 100.0;
 319   cmax = 0.0;
 320   for(i=0; (i<n_x); i++)
 321     for(j=0; (j<max_size); j++) {
 322       if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
 323         cmid = cs_dist[i][j];
 324       cmax = max(cs_dist[i][j],cmax);
 325     }
 326   fprintf(stderr,"cmid: %g, cmax: %g, max_size: %d\n",cmid,cmax,max_size);
 327   cmid = 1;
 328   fp = ffopen(xpm,"w");
 329   write_xpm3(fp,0,"Cluster size distribution","# clusters",timebuf,"Size",
 330              n_x,max_size,t_x,t_y,cs_dist,0,cmid,cmax,
 331              rlo,rmid,rhi,&nlevels);
 332   ffclose(fp);
 333   cmid = 100.0;
 334   cmax = 0.0;
 335   for(i=0; (i<n_x); i++)
 336     for(j=0; (j<max_size); j++) {
 337       cs_dist[i][j] *= (j+1);
 338       if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
 339         cmid = cs_dist[i][j];
 340       cmax = max(cs_dist[i][j],cmax);
 341     }
 342   fprintf(stderr,"cmid: %g, cmax: %g, max_size: %d\n",cmid,cmax,max_size);
 343   fp = ffopen(xpmw,"w");
 344   write_xpm3(fp,0,"Weighted cluster size distribution","Fraction",timebuf,
 345              "Size", n_x,max_size,t_x,t_y,cs_dist,0,cmid,cmax,
 346              rlo,rmid,rhi,&nlevels);
 347   ffclose(fp);
 348
 349   sfree(clust_index);
 350   sfree(clust_size);
 351   sfree(index);
 352 }
 353
 354 int gmx_clustsize(int argc,char *argv[])
 355 {
 356   const char *desc[] = {
 357     "This program computes the size distributions of molecular/atomic clusters in",
 358     "the gas phase. The output is given in the form of a XPM file.",
 359     "The total number of clusters is written to a XVG file.[PAR]",
 360     "When the [TT]-mol[tt] option is given clusters will be made out of",
 361     "molecules rather than atoms, which allows clustering of large molecules.",
 362     "In this case an index file would still contain atom numbers",
 363     "or your calculcation will die with a SEGV.[PAR]",
 364     "When velocities are present in your trajectory, the temperature of",
 365     "the largest cluster will be printed in a separate xvg file assuming",
 366     "that the particles are free to move. If you are using constraints,",
 367     "please correct the temperature. For instance water simulated with SHAKE",
 368     "or SETTLE will yield a temperature that is 1.5 times too low. You can",
 369     "compensate for this with the -ndf option. Remember to take the removal",
 370     "of center of mass motion into account.[PAR]",
 371     "The [TT]-mc[tt] option will produce an index file containing the",
 372     "atom numbers of the largest cluster."
 373   };
 374
 375   static real cutoff   = 0.35;
 376   static int  nskip    = 0;
 377   static int  nlevels  = 20;
 378   static int  ndf      = -1;
 379   static bool bMol     = FALSE;
 380   static bool bPBC     = TRUE;
 381   static rvec rlo      = { 1.0, 1.0, 0.0 };
 382   static rvec rhi      = { 0.0, 0.0, 1.0 };
 383
 384   output_env_t oenv;
 385
 386   t_pargs pa[] = {
 387     { "-cut",      FALSE, etREAL, {&cutoff},
 388       "Largest distance (nm) to be considered in a cluster" },
 389     { "-mol",      FALSE, etBOOL, {&bMol},
 390       "Cluster molecules rather than atoms (needs tpr file)" },
 391     { "-pbc",      FALSE, etBOOL, {&bPBC},
 392       "Use periodic boundary conditions" },
 393     { "-nskip",    FALSE, etINT,  {&nskip},
 394       "Number of frames to skip between writing" },
 395     { "-nlevels",  FALSE, etINT,  {&nlevels},
 396       "Number of levels of grey in xpm output" },
 397     { "-ndf",      FALSE, etINT,  {&ndf},
 398       "Number of degrees of freedom of the entire system for temperature calculation. If not set the number of atoms times three is used." },
 399     { "-rgblo",    FALSE, etRVEC, {rlo},
 400       "RGB values for the color of the lowest occupied cluster size" },
 401     { "-rgbhi",    FALSE, etRVEC, {rhi},
 402       "RGB values for the color of the highest occupied cluster size" }
 403   };
 404 #define NPA asize(pa)
 405   const char *fnNDX,*fnTPR;
 406   bool       bSQ,bRDF;
 407   t_rgb      rgblo,rgbhi;
 408
 409   t_filenm   fnm[] = {
 410     { efTRX, "-f",  NULL,         ffREAD  },
 411     { efTPR, NULL,  NULL,         ffOPTRD },
 412     { efNDX, NULL,  NULL,         ffOPTRD },
 413     { efXPM, "-o", "csize",       ffWRITE },
 414     { efXPM, "-ow","csizew",      ffWRITE },
 415     { efXVG, "-nc","nclust",      ffWRITE },
 416     { efXVG, "-mc","maxclust",    ffWRITE },
 417     { efXVG, "-ac","avclust",     ffWRITE },
 418     { efXVG, "-hc","histo-clust", ffWRITE },
 419     { efXVG, "-temp","temp",     ffOPTWR },
 420     { efNDX, "-mcn", "maxclust", ffOPTWR }
 421   };
 422 #define NFILE asize(fnm)
 423
 424   CopyRight(stderr,argv[0]);
 425   parse_common_args(&argc,argv,
 426                     PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 427                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 428
 429   fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
 430   rgblo.r = rlo[XX],rgblo.g = rlo[YY],rgblo.b = rlo[ZZ];
 431   rgbhi.r = rhi[XX],rgbhi.g = rhi[YY],rgbhi.b = rhi[ZZ];
 432
 433   fnTPR = ftp2fn_null(efTPR,NFILE,fnm);
 434   if (bMol && !fnTPR)
 435     gmx_fatal(FARGS,"You need a tpr file for the -mol option");
 436
 437   clust_size(fnNDX,ftp2fn(efTRX,NFILE,fnm),opt2fn("-o",NFILE,fnm),
 438              opt2fn("-ow",NFILE,fnm),
 439              opt2fn("-nc",NFILE,fnm),opt2fn("-ac",NFILE,fnm),
 440              opt2fn("-mc",NFILE,fnm),opt2fn("-hc",NFILE,fnm),
 441              opt2fn("-temp",NFILE,fnm),opt2fn("-mcn",NFILE,fnm),
 442              bMol,bPBC,fnTPR,
 443              cutoff,nskip,nlevels,rgblo,rgbhi,ndf,oenv);
 444
 445   thanx(stderr);
 446
 447   return 0;
 448 }