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