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