src/gromacs/gmxana/gmx_clustsize.c

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