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, "negative cluster size %d for element %d",
 270                                                   clust_size[cj], cj);
 271                                     }
 272                                     clust_size[cj]--;
 273                                     clust_index[k] = ci;
 274                                     clust_size[ci]++;
 275                                 }
 276                             }
 277                         }
 278                     }
 279                 }
 280             }
 281             n_x++;
 282             srenew(t_x, n_x);
 283             if (fr.bTime)
 284             {
 285                 frameTime = fr.time;
 286             }
 287             else if (fr.bStep)
 288             {
 289                 frameTime = fr.step;
 290             }
 291             else
 292             {
 293                 frameTime = ++frameCounter;
 294             }
 295             t_x[n_x - 1] = frameTime * tf;
 296             srenew(cs_dist, n_x);
 297             snew(cs_dist[n_x - 1], nindex);
 298             nclust = 0;
 299             cav    = 0;
 300             nav    = 0;
 301             for (i = 0; (i < nindex); i++)
 302             {
 303                 ci = clust_size[i];
 304                 if (ci > max_clust_size)
 305                 {
 306                     max_clust_size = ci;
 307                     max_clust_ind  = i;
 308                 }
 309                 if (ci > 0)
 310                 {
 311                     nclust++;
 312                     cs_dist[n_x - 1][ci - 1] += 1.0;
 313                     max_size = std::max(max_size, ci);
 314                     if (ci > 1)
 315                     {
 316                         cav += ci;
 317                         nav++;
 318                     }
 319                 }
 320             }
 321             fprintf(fp, "%14.6e  %10d\n", frameTime, nclust);
 322             if (nav > 0)
 323             {
 324                 fprintf(gp, "%14.6e  %10.3f\n", frameTime, cav / nav);
 325             }
 326             fprintf(hp, "%14.6e  %10d\n", frameTime, max_clust_size);
 327         }
 328         /* Analyse velocities, if present */
 329         if (fr.bV)
 330         {
 331             if (!tpr)
 332             {
 333                 if (bTPRwarn)
 334                 {
 335                     printf("You need a [REF].tpr[ref] file to analyse temperatures\n");
 336                     bTPRwarn = FALSE;
 337                 }
 338             }
 339             else
 340             {
 341                 v = fr.v;
 342                 /* Loop over clusters and for each cluster compute 1/2 m v^2 */
 343                 if (max_clust_ind >= 0)
 344                 {
 345                     ekin = 0;
 346                     for (i = 0; (i < nindex); i++)
 347                     {
 348                         if (clust_index[i] == max_clust_ind)
 349                         {
 350                             ai     = index[i];
 351                             real m = mtopGetAtomMass(mtop, ai, &molb);
 352                             ekin += 0.5 * m * iprod(v[ai], v[ai]);
 353                         }
 354                     }
 355                     temp = (ekin * 2.0) / (3.0 * tfac * max_clust_size * BOLTZ);
 356                     fprintf(tp, "%10.3f  %10.3f\n", frameTime, temp);
 357                 }
 358             }
 359         }
 360         nframe++;
 361     } while (read_next_frame(oenv, status, &fr));
 362     close_trx(status);
 363     done_frame(&fr);
 364     xvgrclose(fp);
 365     xvgrclose(gp);
 366     xvgrclose(hp);
 367     xvgrclose(tp);
 368
 369     if (max_clust_ind >= 0)
 370     {
 371         fp = gmx_ffopen(mcn, "w");
 372         fprintf(fp, "[ max_clust ]\n");
 373         for (i = 0; (i < nindex); i++)
 374         {
 375             if (clust_index[i] == max_clust_ind)
 376             {
 377                 if (bMol)
 378                 {
 379                     GMX_RELEASE_ASSERT(mols.numBlocks() > 0,
 380                                        "Cannot access index[] from empty mols");
 381                     for (int j : mols.block(i))
 382                     {
 383                         fprintf(fp, "%d\n", j + 1);
 384                     }
 385                 }
 386                 else
 387                 {
 388                     fprintf(fp, "%d\n", index[i] + 1);
 389                 }
 390             }
 391         }
 392         gmx_ffclose(fp);
 393     }
 394
 395     /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
 396     fp     = xvgropen(histo, "Cluster size distribution", "Cluster size", "()", oenv);
 397     nhisto = 0;
 398     fprintf(fp, "%5d  %8.3f\n", 0, 0.0);
 399     for (j = 0; (j < max_size); j++)
 400     {
 401         real nelem = 0;
 402         for (i = 0; (i < n_x); i++)
 403         {
 404             nelem += cs_dist[i][j];
 405         }
 406         fprintf(fp, "%5d  %8.3f\n", j + 1, nelem / n_x);
 407         nhisto += static_cast<int>((j + 1) * nelem / n_x);
 408     }
 409     fprintf(fp, "%5d  %8.3f\n", j + 1, 0.0);
 410     xvgrclose(fp);
 411
 412     fprintf(stderr, "Total number of atoms in clusters =  %d\n", nhisto);
 413
 414     /* Look for the smallest entry that is not zero
 415      * This will make that zero is white, and not zero is coloured.
 416      */
 417     cmid = 100.0;
 418     cmax = 0.0;
 419     for (i = 0; (i < n_x); i++)
 420     {
 421         for (j = 0; (j < max_size); j++)
 422         {
 423             if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
 424             {
 425                 cmid = cs_dist[i][j];
 426             }
 427             cmax = std::max(cs_dist[i][j], cmax);
 428         }
 429     }
 430     fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
 431     cmid = 1;
 432     fp   = gmx_ffopen(xpm, "w");
 433     write_xpm3(fp, 0, "Cluster size distribution", "# clusters", timeLabel, "Size", n_x, max_size,
 434                t_x, t_y, cs_dist, 0, cmid, cmax, rlo, rmid, rhi, &nlevels);
 435     gmx_ffclose(fp);
 436     cmid = 100.0;
 437     cmax = 0.0;
 438     for (i = 0; (i < n_x); i++)
 439     {
 440         for (j = 0; (j < max_size); j++)
 441         {
 442             cs_dist[i][j] *= (j + 1);
 443             if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
 444             {
 445                 cmid = cs_dist[i][j];
 446             }
 447             cmax = std::max(cs_dist[i][j], cmax);
 448         }
 449     }
 450     fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
 451     fp = gmx_ffopen(xpmw, "w");
 452     write_xpm3(fp, 0, "Weighted cluster size distribution", "Fraction", timeLabel, "Size", n_x,
 453                max_size, t_x, t_y, cs_dist, 0, cmid, cmax, rlo, rmid, rhi, &nlevels);
 454     gmx_ffclose(fp);
 455     delete mtop;
 456     sfree(t_x);
 457     sfree(t_y);
 458     for (i = 0; (i < n_x); i++)
 459     {
 460         sfree(cs_dist[i]);
 461     }
 462     sfree(cs_dist);
 463     sfree(clust_index);
 464     sfree(clust_size);
 465     sfree(index);
 466 }
 467
 468 int gmx_clustsize(int argc, char* argv[])
 469 {
 470     const char* desc[] = {
 471         "[THISMODULE] computes the size distributions of molecular/atomic clusters in",
 472         "the gas phase. The output is given in the form of an [REF].xpm[ref] file.",
 473         "The total number of clusters is written to an [REF].xvg[ref] file.[PAR]",
 474         "When the [TT]-mol[tt] option is given clusters will be made out of",
 475         "molecules rather than atoms, which allows clustering of large molecules.",
 476         "In this case an index file would still contain atom numbers",
 477         "or your calculation will die with a SEGV.[PAR]",
 478         "When velocities are present in your trajectory, the temperature of",
 479         "the largest cluster will be printed in a separate [REF].xvg[ref] file assuming",
 480         "that the particles are free to move. If you are using constraints,",
 481         "please correct the temperature. For instance water simulated with SHAKE",
 482         "or SETTLE will yield a temperature that is 1.5 times too low. You can",
 483         "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal",
 484         "of center of mass motion into account.[PAR]",
 485         "The [TT]-mc[tt] option will produce an index file containing the",
 486         "atom numbers of the largest cluster."
 487     };
 488
 489     real     cutoff  = 0.35;
 490     int      nskip   = 0;
 491     int      nlevels = 20;
 492     int      ndf     = -1;
 493     gmx_bool bMol    = FALSE;
 494     gmx_bool bPBC    = TRUE;
 495     rvec     rlo     = { 1.0, 1.0, 0.0 };
 496     rvec     rhi     = { 0.0, 0.0, 1.0 };
 497
 498     gmx_output_env_t* oenv;
 499
 500     t_pargs pa[] = {
 501         { "-cut",
 502           FALSE,
 503           etREAL,
 504           { &cutoff },
 505           "Largest distance (nm) to be considered in a cluster" },
 506         { "-mol",
 507           FALSE,
 508           etBOOL,
 509           { &bMol },
 510           "Cluster molecules rather than atoms (needs [REF].tpr[ref] file)" },
 511         { "-pbc", FALSE, etBOOL, { &bPBC }, "Use periodic boundary conditions" },
 512         { "-nskip", FALSE, etINT, { &nskip }, "Number of frames to skip between writing" },
 513         { "-nlevels",
 514           FALSE,
 515           etINT,
 516           { &nlevels },
 517           "Number of levels of grey in [REF].xpm[ref] output" },
 518         { "-ndf",
 519           FALSE,
 520           etINT,
 521           { &ndf },
 522           "Number of degrees of freedom of the entire system for temperature calculation. "
 523           "If not set, the number of atoms times three is used." },
 524         { "-rgblo",
 525           FALSE,
 526           etRVEC,
 527           { rlo },
 528           "RGB values for the color of the lowest occupied cluster size" },
 529         { "-rgbhi",
 530           FALSE,
 531           etRVEC,
 532           { rhi },
 533           "RGB values for the color of the highest occupied cluster size" }
 534     };
 535 #define NPA asize(pa)
 536     const char *fnNDX, *fnTPR;
 537     t_rgb       rgblo, rgbhi;
 538
 539     t_filenm fnm[] = {
 540         { efTRX, "-f", nullptr, ffREAD },         { efTPR, nullptr, nullptr, ffOPTRD },
 541         { efNDX, nullptr, nullptr, ffOPTRD },     { efXPM, "-o", "csize", ffWRITE },
 542         { efXPM, "-ow", "csizew", ffWRITE },      { efXVG, "-nc", "nclust", ffWRITE },
 543         { efXVG, "-mc", "maxclust", ffWRITE },    { efXVG, "-ac", "avclust", ffWRITE },
 544         { efXVG, "-hc", "histo-clust", ffWRITE }, { efXVG, "-temp", "temp", ffOPTWR },
 545         { efNDX, "-mcn", "maxclust", ffOPTWR }
 546     };
 547 #define NFILE asize(fnm)
 548
 549     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm,
 550                            NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 551     {
 552         return 0;
 553     }
 554
 555     fnNDX   = ftp2fn_null(efNDX, NFILE, fnm);
 556     rgblo.r = rlo[XX];
 557     rgblo.g = rlo[YY];
 558     rgblo.b = rlo[ZZ];
 559     rgbhi.r = rhi[XX];
 560     rgbhi.g = rhi[YY];
 561     rgbhi.b = rhi[ZZ];
 562
 563     fnTPR = ftp2fn_null(efTPR, NFILE, fnm);
 564     if (bMol && !fnTPR)
 565     {
 566         gmx_fatal(FARGS, "You need a tpr file for the -mol option");
 567     }
 568
 569     clust_size(fnNDX, ftp2fn(efTRX, NFILE, fnm), opt2fn("-o", NFILE, fnm), opt2fn("-ow", NFILE, fnm),
 570                opt2fn("-nc", NFILE, fnm), opt2fn("-ac", NFILE, fnm), opt2fn("-mc", NFILE, fnm),
 571                opt2fn("-hc", NFILE, fnm), opt2fn("-temp", NFILE, fnm), opt2fn("-mcn", NFILE, fnm),
 572                bMol, bPBC, fnTPR, cutoff, nskip, nlevels, rgblo, rgbhi, ndf, oenv);
 573
 574     output_env_done(oenv);
 575
 576     return 0;
 577 }