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