src/gromacs/gmxana/gmx_clustsize.cpp

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