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