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