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