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