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