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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source 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 #include "gmxpre.h"
  39
  40 #include <cmath>
  41
  42 #include <algorithm>
  43
  44 #include "gromacs/commandline/filenm.h"
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/fileio/matio.h"
  47 #include "gromacs/fileio/tpxio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/gmx_ana.h"
  51 #include "gromacs/gmxana/gstat.h"
  52 #include "gromacs/gmxlib/nrnb.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/pbcutil/pbc.h"
  56 #include "gromacs/topology/index.h"
  57 #include "gromacs/topology/mtop_lookup.h"
  58 #include "gromacs/topology/mtop_util.h"
  59 #include "gromacs/topology/topology.h"
  60 #include "gromacs/trajectory/trajectoryframe.h"
  61 #include "gromacs/utility/arraysize.h"
  62 #include "gromacs/utility/cstringutil.h"
  63 #include "gromacs/utility/fatalerror.h"
  64 #include "gromacs/utility/futil.h"
  65 #include "gromacs/utility/gmxassert.h"
  66 #include "gromacs/utility/smalloc.h"
  67
  68 static void clust_size(const char*             ndx,
  69                        const char*             trx,
  70                        const char*             xpm,
  71                        const char*             xpmw,
  72                        const char*             ncl,
  73                        const char*             acl,
  74                        const char*             mcl,
  75                        const char*             histo,
  76                        const char*             tempf,
  77                        const char*             mcn,
  78                        gmx_bool                bMol,
  79                        gmx_bool                bPBC,
  80                        const char*             tpr,
  81                        real                    cut,
  82                        int                     nskip,
  83                        int                     nlevels,
  84                        t_rgb                   rmid,
  85                        t_rgb                   rhi,
  86                        int                     ndf,
  87                        const gmx_output_env_t* oenv)
  88 {
  89     FILE *       fp, *gp, *hp, *tp;
  90     int*         index = nullptr;
  91     int          nindex, natoms;
  92     t_trxstatus* status;
  93     rvec *       x = nullptr, *v = nullptr, dx;
  94     t_pbc        pbc;
  95     gmx_bool     bSame, bTPRwarn = TRUE;
  96     /* Topology stuff */
  97     t_trxframe    fr;
  98     TpxFileHeader tpxh;
  99     gmx_mtop_t    mtop;
 100     PbcType       pbcType = PbcType::Unset;
 101     int           ii, jj;
 102     real          temp, tfac;
 103     /* Cluster size distribution (matrix) */
 104     real** cs_dist = nullptr;
 105     real   tf, dx2, cut2, *t_x = nullptr, *t_y, cmid, cmax, cav, ekin;
 106     int    i, j, k, ai, aj, ci, cj, nframe, nclust, n_x, max_size = 0;
 107     int *  clust_index, *clust_size, max_clust_size, max_clust_ind, nav, nhisto;
 108     t_rgb  rlo          = { 1.0, 1.0, 1.0 };
 109     int    frameCounter = 0;
 110     real   frameTime;
 111
 112     clear_trxframe(&fr, TRUE);
 113     auto timeLabel = output_env_get_time_label(oenv);
 114     tf             = output_env_get_time_factor(oenv);
 115     fp             = xvgropen(ncl, "Number of clusters", timeLabel, "N", oenv);
 116     gp             = xvgropen(acl, "Average cluster size", timeLabel, "#molecules", oenv);
 117     hp             = xvgropen(mcl, "Max cluster size", timeLabel, "#molecules", oenv);
 118     tp             = xvgropen(tempf, "Temperature of largest cluster", timeLabel, "T (K)", oenv);
 119
 120     if (!read_first_frame(oenv, &status, trx, &fr, TRX_NEED_X | TRX_READ_V))
 121     {
 122         gmx_file(trx);
 123     }
 124
 125     natoms = fr.natoms;
 126     x      = fr.x;
 127
 128     if (tpr)
 129     {
 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         pbcType = 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(tpr, "Cannot access topology without having read it from TPR");
 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, pbcType, 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,
 269                                                   "negative cluster size %d for element %d",
 270                                                   clust_size[cj],
 271                                                   cj);
 272                                     }
 273                                     clust_size[cj]--;
 274                                     clust_index[k] = ci;
 275                                     clust_size[ci]++;
 276                                 }
 277                             }
 278                         }
 279                     }
 280                 }
 281             }
 282             n_x++;
 283             srenew(t_x, n_x);
 284             if (fr.bTime)
 285             {
 286                 frameTime = fr.time;
 287             }
 288             else if (fr.bStep)
 289             {
 290                 frameTime = fr.step;
 291             }
 292             else
 293             {
 294                 frameTime = ++frameCounter;
 295             }
 296             t_x[n_x - 1] = frameTime * tf;
 297             srenew(cs_dist, n_x);
 298             snew(cs_dist[n_x - 1], nindex);
 299             nclust = 0;
 300             cav    = 0;
 301             nav    = 0;
 302             for (i = 0; (i < nindex); i++)
 303             {
 304                 ci = clust_size[i];
 305                 if (ci > max_clust_size)
 306                 {
 307                     max_clust_size = ci;
 308                     max_clust_ind  = i;
 309                 }
 310                 if (ci > 0)
 311                 {
 312                     nclust++;
 313                     cs_dist[n_x - 1][ci - 1] += 1.0;
 314                     max_size = std::max(max_size, ci);
 315                     if (ci > 1)
 316                     {
 317                         cav += ci;
 318                         nav++;
 319                     }
 320                 }
 321             }
 322             fprintf(fp, "%14.6e  %10d\n", frameTime, nclust);
 323             if (nav > 0)
 324             {
 325                 fprintf(gp, "%14.6e  %10.3f\n", frameTime, cav / nav);
 326             }
 327             fprintf(hp, "%14.6e  %10d\n", frameTime, max_clust_size);
 328         }
 329         /* Analyse velocities, if present */
 330         if (fr.bV)
 331         {
 332             if (!tpr)
 333             {
 334                 if (bTPRwarn)
 335                 {
 336                     printf("You need a [REF].tpr[ref] file to analyse temperatures\n");
 337                     bTPRwarn = FALSE;
 338                 }
 339             }
 340             else
 341             {
 342                 v = fr.v;
 343                 /* Loop over clusters and for each cluster compute 1/2 m v^2 */
 344                 if (max_clust_ind >= 0)
 345                 {
 346                     ekin = 0;
 347                     for (i = 0; (i < nindex); i++)
 348                     {
 349                         if (clust_index[i] == max_clust_ind)
 350                         {
 351                             ai     = index[i];
 352                             real m = mtopGetAtomMass(mtop, ai, &molb);
 353                             ekin += 0.5 * m * iprod(v[ai], v[ai]);
 354                         }
 355                     }
 356                     temp = (ekin * 2.0) / (3.0 * tfac * max_clust_size * gmx::c_boltz);
 357                     fprintf(tp, "%10.3f  %10.3f\n", frameTime, temp);
 358                 }
 359             }
 360         }
 361         nframe++;
 362     } while (read_next_frame(oenv, status, &fr));
 363     close_trx(status);
 364     done_frame(&fr);
 365     xvgrclose(fp);
 366     xvgrclose(gp);
 367     xvgrclose(hp);
 368     xvgrclose(tp);
 369
 370     if (max_clust_ind >= 0)
 371     {
 372         fp = gmx_ffopen(mcn, "w");
 373         fprintf(fp, "[ max_clust ]\n");
 374         for (i = 0; (i < nindex); i++)
 375         {
 376             if (clust_index[i] == max_clust_ind)
 377             {
 378                 if (bMol)
 379                 {
 380                     GMX_RELEASE_ASSERT(mols.numBlocks() > 0,
 381                                        "Cannot access index[] from empty mols");
 382                     for (int j : mols.block(i))
 383                     {
 384                         fprintf(fp, "%d\n", j + 1);
 385                     }
 386                 }
 387                 else
 388                 {
 389                     fprintf(fp, "%d\n", index[i] + 1);
 390                 }
 391             }
 392         }
 393         gmx_ffclose(fp);
 394     }
 395
 396     /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
 397     fp     = xvgropen(histo, "Cluster size distribution", "Cluster size", "()", oenv);
 398     nhisto = 0;
 399     fprintf(fp, "%5d  %8.3f\n", 0, 0.0);
 400     for (j = 0; (j < max_size); j++)
 401     {
 402         real nelem = 0;
 403         for (i = 0; (i < n_x); i++)
 404         {
 405             nelem += cs_dist[i][j];
 406         }
 407         fprintf(fp, "%5d  %8.3f\n", j + 1, nelem / n_x);
 408         nhisto += static_cast<int>((j + 1) * nelem / n_x);
 409     }
 410     fprintf(fp, "%5d  %8.3f\n", j + 1, 0.0);
 411     xvgrclose(fp);
 412
 413     fprintf(stderr, "Total number of atoms in clusters =  %d\n", nhisto);
 414
 415     /* Look for the smallest entry that is not zero
 416      * This will make that zero is white, and not zero is coloured.
 417      */
 418     cmid = 100.0;
 419     cmax = 0.0;
 420     for (i = 0; (i < n_x); i++)
 421     {
 422         for (j = 0; (j < max_size); j++)
 423         {
 424             if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
 425             {
 426                 cmid = cs_dist[i][j];
 427             }
 428             cmax = std::max(cs_dist[i][j], cmax);
 429         }
 430     }
 431     fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
 432     cmid = 1;
 433     fp   = gmx_ffopen(xpm, "w");
 434     write_xpm3(fp,
 435                0,
 436                "Cluster size distribution",
 437                "# clusters",
 438                timeLabel,
 439                "Size",
 440                n_x,
 441                max_size,
 442                t_x,
 443                t_y,
 444                cs_dist,
 445                0,
 446                cmid,
 447                cmax,
 448                rlo,
 449                rmid,
 450                rhi,
 451                &nlevels);
 452     gmx_ffclose(fp);
 453     cmid = 100.0;
 454     cmax = 0.0;
 455     for (i = 0; (i < n_x); i++)
 456     {
 457         for (j = 0; (j < max_size); j++)
 458         {
 459             cs_dist[i][j] *= (j + 1);
 460             if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
 461             {
 462                 cmid = cs_dist[i][j];
 463             }
 464             cmax = std::max(cs_dist[i][j], cmax);
 465         }
 466     }
 467     fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
 468     fp = gmx_ffopen(xpmw, "w");
 469     write_xpm3(fp,
 470                0,
 471                "Weighted cluster size distribution",
 472                "Fraction",
 473                timeLabel,
 474                "Size",
 475                n_x,
 476                max_size,
 477                t_x,
 478                t_y,
 479                cs_dist,
 480                0,
 481                cmid,
 482                cmax,
 483                rlo,
 484                rmid,
 485                rhi,
 486                &nlevels);
 487     gmx_ffclose(fp);
 488     sfree(t_x);
 489     sfree(t_y);
 490     for (i = 0; (i < n_x); i++)
 491     {
 492         sfree(cs_dist[i]);
 493     }
 494     sfree(cs_dist);
 495     sfree(clust_index);
 496     sfree(clust_size);
 497     sfree(index);
 498 }
 499
 500 int gmx_clustsize(int argc, char* argv[])
 501 {
 502     const char* desc[] = {
 503         "[THISMODULE] computes the size distributions of molecular/atomic clusters in",
 504         "the gas phase. The output is given in the form of an [REF].xpm[ref] file.",
 505         "The total number of clusters is written to an [REF].xvg[ref] file.[PAR]",
 506         "When the [TT]-mol[tt] option is given clusters will be made out of",
 507         "molecules rather than atoms, which allows clustering of large molecules.",
 508         "In this case an index file would still contain atom numbers",
 509         "or your calculation will die with a SEGV.[PAR]",
 510         "When velocities are present in your trajectory, the temperature of",
 511         "the largest cluster will be printed in a separate [REF].xvg[ref] file assuming",
 512         "that the particles are free to move. If you are using constraints,",
 513         "please correct the temperature. For instance water simulated with SHAKE",
 514         "or SETTLE will yield a temperature that is 1.5 times too low. You can",
 515         "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal",
 516         "of center of mass motion into account.[PAR]",
 517         "The [TT]-mc[tt] option will produce an index file containing the",
 518         "atom numbers of the largest cluster."
 519     };
 520
 521     real     cutoff  = 0.35;
 522     int      nskip   = 0;
 523     int      nlevels = 20;
 524     int      ndf     = -1;
 525     gmx_bool bMol    = FALSE;
 526     gmx_bool bPBC    = TRUE;
 527     rvec     rlo     = { 1.0, 1.0, 0.0 };
 528     rvec     rhi     = { 0.0, 0.0, 1.0 };
 529
 530     gmx_output_env_t* oenv;
 531
 532     t_pargs pa[] = {
 533         { "-cut",
 534           FALSE,
 535           etREAL,
 536           { &cutoff },
 537           "Largest distance (nm) to be considered in a cluster" },
 538         { "-mol",
 539           FALSE,
 540           etBOOL,
 541           { &bMol },
 542           "Cluster molecules rather than atoms (needs [REF].tpr[ref] file)" },
 543         { "-pbc", FALSE, etBOOL, { &bPBC }, "Use periodic boundary conditions" },
 544         { "-nskip", FALSE, etINT, { &nskip }, "Number of frames to skip between writing" },
 545         { "-nlevels",
 546           FALSE,
 547           etINT,
 548           { &nlevels },
 549           "Number of levels of grey in [REF].xpm[ref] output" },
 550         { "-ndf",
 551           FALSE,
 552           etINT,
 553           { &ndf },
 554           "Number of degrees of freedom of the entire system for temperature calculation. "
 555           "If not set, the number of atoms times three is used." },
 556         { "-rgblo",
 557           FALSE,
 558           etRVEC,
 559           { rlo },
 560           "RGB values for the color of the lowest occupied cluster size" },
 561         { "-rgbhi",
 562           FALSE,
 563           etRVEC,
 564           { rhi },
 565           "RGB values for the color of the highest occupied cluster size" }
 566     };
 567 #define NPA asize(pa)
 568     const char *fnNDX, *fnTPR;
 569     t_rgb       rgblo, rgbhi;
 570
 571     t_filenm fnm[] = {
 572         { efTRX, "-f", nullptr, ffREAD },         { efTPR, nullptr, nullptr, ffOPTRD },
 573         { efNDX, nullptr, nullptr, ffOPTRD },     { efXPM, "-o", "csize", ffWRITE },
 574         { efXPM, "-ow", "csizew", ffWRITE },      { efXVG, "-nc", "nclust", ffWRITE },
 575         { efXVG, "-mc", "maxclust", ffWRITE },    { efXVG, "-ac", "avclust", ffWRITE },
 576         { efXVG, "-hc", "histo-clust", ffWRITE }, { efXVG, "-temp", "temp", ffOPTWR },
 577         { efNDX, "-mcn", "maxclust", ffOPTWR }
 578     };
 579 #define NFILE asize(fnm)
 580
 581     if (!parse_common_args(
 582                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 583     {
 584         return 0;
 585     }
 586
 587     fnNDX   = ftp2fn_null(efNDX, NFILE, fnm);
 588     rgblo.r = rlo[XX];
 589     rgblo.g = rlo[YY];
 590     rgblo.b = rlo[ZZ];
 591     rgbhi.r = rhi[XX];
 592     rgbhi.g = rhi[YY];
 593     rgbhi.b = rhi[ZZ];
 594
 595     fnTPR = ftp2fn_null(efTPR, NFILE, fnm);
 596     if (bMol && !fnTPR)
 597     {
 598         gmx_fatal(FARGS, "You need a tpr file for the -mol option");
 599     }
 600
 601     clust_size(fnNDX,
 602                ftp2fn(efTRX, NFILE, fnm),
 603                opt2fn("-o", NFILE, fnm),
 604                opt2fn("-ow", NFILE, fnm),
 605                opt2fn("-nc", NFILE, fnm),
 606                opt2fn("-ac", NFILE, fnm),
 607                opt2fn("-mc", NFILE, fnm),
 608                opt2fn("-hc", NFILE, fnm),
 609                opt2fn("-temp", NFILE, fnm),
 610                opt2fn("-mcn", NFILE, fnm),
 611                bMol,
 612                bPBC,
 613                fnTPR,
 614                cutoff,
 615                nskip,
 616                nlevels,
 617                rgblo,
 618                rgbhi,
 619                ndf,
 620                oenv);
 621
 622     output_env_done(oenv);
 623
 624     return 0;
 625 }