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