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