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