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