src/tools/gmx_clustsize.c

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