3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2007, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Groningen Machine for Chemical Simulation
63 #include "mtop_util.h"
67 static void clust_size(const char *ndx, const char *trx, const char *xpm,
68 const char *xpmw, const char *ncl, const char *acl,
69 const char *mcl, const char *histo, const char *tempf,
70 const char *mcn, gmx_bool bMol, gmx_bool bPBC, const char *tpr,
71 real cut, int nskip, int nlevels,
72 t_rgb rmid, t_rgb rhi, int ndf,
73 const output_env_t oenv)
75 FILE *fp, *gp, *hp, *tp;
76 atom_id *index = NULL;
79 rvec *x = NULL, *v = NULL, dx;
83 gmx_bool bSame, bTPRwarn = TRUE;
87 gmx_mtop_t *mtop = NULL;
90 gmx_mtop_atomlookup_t alook;
92 int version, generation, ii, jj, nsame;
94 /* Cluster size distribution (matrix) */
95 real **cs_dist = NULL;
96 real tf, dx2, cut2, *t_x = NULL, *t_y, cmid, cmax, cav, ekin;
97 int i, j, k, ai, aj, ak, ci, cj, nframe, nclust, n_x, n_y, max_size = 0;
98 int *clust_index, *clust_size, max_clust_size, max_clust_ind, nav, nhisto;
99 t_rgb rlo = { 1.0, 1.0, 1.0 };
101 clear_trxframe(&fr, TRUE);
102 sprintf(timebuf, "Time (%s)", output_env_get_time_unit(oenv));
103 tf = output_env_get_time_factor(oenv);
104 fp = xvgropen(ncl, "Number of clusters", timebuf, "N", oenv);
105 gp = xvgropen(acl, "Average cluster size", timebuf, "#molecules", oenv);
106 hp = xvgropen(mcl, "Max cluster size", timebuf, "#molecules", oenv);
107 tp = xvgropen(tempf, "Temperature of largest cluster", timebuf, "T (K)",
110 if (!read_first_frame(oenv, &status, trx, &fr, TRX_NEED_X | TRX_READ_V))
121 read_tpxheader(tpr, &tpxh, TRUE, &version, &generation);
122 if (tpxh.natoms != natoms)
124 gmx_fatal(FARGS, "tpr (%d atoms) and xtc (%d atoms) do not match!",
125 tpxh.natoms, natoms);
127 ePBC = read_tpx(tpr, NULL, NULL, &natoms, NULL, NULL, NULL, mtop);
135 tfac = ndf/(3.0*natoms);
142 printf("Using molecules rather than atoms. Not reading index file %s\n",
145 mols = &(mtop->mols);
147 /* Make dummy index */
150 for (i = 0; (i < nindex); i++)
154 gname = strdup("mols");
158 rd_index(ndx, 1, &nindex, &index, &gname);
161 alook = gmx_mtop_atomlookup_init(mtop);
163 snew(clust_index, nindex);
164 snew(clust_size, nindex);
169 for (i = 0; (i < nindex); i++)
177 if ((nskip == 0) || ((nskip > 0) && ((nframe % nskip) == 0)))
181 set_pbc(&pbc, ePBC, fr.box);
186 /* Put all atoms/molecules in their own cluster, with size 1 */
187 for (i = 0; (i < nindex); i++)
189 /* Cluster index is indexed with atom index number */
191 /* Cluster size is indexed with cluster number */
195 /* Loop over atoms */
196 for (i = 0; (i < nindex); i++)
201 /* Loop over atoms (only half a matrix) */
202 for (j = i+1; (j < nindex); j++)
206 /* If they are not in the same cluster already */
211 /* Compute distance */
215 for (ii = mols->index[ai]; !bSame && (ii < mols->index[ai+1]); ii++)
217 for (jj = mols->index[aj]; !bSame && (jj < mols->index[aj+1]); jj++)
221 pbc_dx(&pbc, x[ii], x[jj], dx);
225 rvec_sub(x[ii], x[jj], dx);
228 bSame = (dx2 < cut2);
236 pbc_dx(&pbc, x[ai], x[aj], dx);
240 rvec_sub(x[ai], x[aj], dx);
243 bSame = (dx2 < cut2);
245 /* If distance less than cut-off */
248 /* Merge clusters: check for all atoms whether they are in
249 * cluster cj and if so, put them in ci
251 for (k = 0; (k < nindex); k++)
253 if (clust_index[k] == cj)
255 if (clust_size[cj] <= 0)
257 gmx_fatal(FARGS, "negative cluster size %d for element %d",
271 t_x[n_x-1] = fr.time*tf;
272 srenew(cs_dist, n_x);
273 snew(cs_dist[n_x-1], nindex);
277 for (i = 0; (i < nindex); i++)
280 if (ci > max_clust_size)
288 cs_dist[n_x-1][ci-1] += 1.0;
289 max_size = max(max_size, ci);
297 fprintf(fp, "%14.6e %10d\n", fr.time, nclust);
300 fprintf(gp, "%14.6e %10.3f\n", fr.time, cav/nav);
302 fprintf(hp, "%14.6e %10d\n", fr.time, max_clust_size);
304 /* Analyse velocities, if present */
311 printf("You need a [TT].tpr[tt] file to analyse temperatures\n");
318 /* Loop over clusters and for each cluster compute 1/2 m v^2 */
319 if (max_clust_ind >= 0)
322 for (i = 0; (i < nindex); i++)
324 if (clust_index[i] == max_clust_ind)
327 gmx_mtop_atomnr_to_atom(alook, ai, &atom);
328 ekin += 0.5*atom->m*iprod(v[ai], v[ai]);
331 temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
332 fprintf(tp, "%10.3f %10.3f\n", fr.time, temp);
338 while (read_next_frame(oenv, status, &fr));
345 gmx_mtop_atomlookup_destroy(alook);
347 if (max_clust_ind >= 0)
349 fp = ffopen(mcn, "w");
350 fprintf(fp, "[ max_clust ]\n");
351 for (i = 0; (i < nindex); i++)
353 if (clust_index[i] == max_clust_ind)
357 for (j = mols->index[i]; (j < mols->index[i+1]); j++)
359 fprintf(fp, "%d\n", j+1);
364 fprintf(fp, "%d\n", index[i]+1);
371 /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
372 fp = xvgropen(histo, "Cluster size distribution", "Cluster size", "()", oenv);
374 fprintf(fp, "%5d %8.3f\n", 0, 0.0);
375 for (j = 0; (j < max_size); j++)
378 for (i = 0; (i < n_x); i++)
380 nelem += cs_dist[i][j];
382 fprintf(fp, "%5d %8.3f\n", j+1, nelem/n_x);
383 nhisto += (int)((j+1)*nelem/n_x);
385 fprintf(fp, "%5d %8.3f\n", j+1, 0.0);
388 fprintf(stderr, "Total number of atoms in clusters = %d\n", nhisto);
390 /* Look for the smallest entry that is not zero
391 * This will make that zero is white, and not zero is coloured.
395 for (i = 0; (i < n_x); i++)
397 for (j = 0; (j < max_size); j++)
399 if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
401 cmid = cs_dist[i][j];
403 cmax = max(cs_dist[i][j], cmax);
406 fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
408 fp = ffopen(xpm, "w");
409 write_xpm3(fp, 0, "Cluster size distribution", "# clusters", timebuf, "Size",
410 n_x, max_size, t_x, t_y, cs_dist, 0, cmid, cmax,
411 rlo, rmid, rhi, &nlevels);
415 for (i = 0; (i < n_x); i++)
417 for (j = 0; (j < max_size); j++)
419 cs_dist[i][j] *= (j+1);
420 if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
422 cmid = cs_dist[i][j];
424 cmax = max(cs_dist[i][j], cmax);
427 fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
428 fp = ffopen(xpmw, "w");
429 write_xpm3(fp, 0, "Weighted cluster size distribution", "Fraction", timebuf,
430 "Size", n_x, max_size, t_x, t_y, cs_dist, 0, cmid, cmax,
431 rlo, rmid, rhi, &nlevels);
439 int gmx_clustsize(int argc, char *argv[])
441 const char *desc[] = {
442 "This program computes the size distributions of molecular/atomic clusters in",
443 "the gas phase. The output is given in the form of an [TT].xpm[tt] file.",
444 "The total number of clusters is written to an [TT].xvg[tt] file.[PAR]",
445 "When the [TT]-mol[tt] option is given clusters will be made out of",
446 "molecules rather than atoms, which allows clustering of large molecules.",
447 "In this case an index file would still contain atom numbers",
448 "or your calculation will die with a SEGV.[PAR]",
449 "When velocities are present in your trajectory, the temperature of",
450 "the largest cluster will be printed in a separate [TT].xvg[tt] file assuming",
451 "that the particles are free to move. If you are using constraints,",
452 "please correct the temperature. For instance water simulated with SHAKE",
453 "or SETTLE will yield a temperature that is 1.5 times too low. You can",
454 "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal",
455 "of center of mass motion into account.[PAR]",
456 "The [TT]-mc[tt] option will produce an index file containing the",
457 "atom numbers of the largest cluster."
460 static real cutoff = 0.35;
461 static int nskip = 0;
462 static int nlevels = 20;
464 static gmx_bool bMol = FALSE;
465 static gmx_bool bPBC = TRUE;
466 static rvec rlo = { 1.0, 1.0, 0.0 };
467 static rvec rhi = { 0.0, 0.0, 1.0 };
472 { "-cut", FALSE, etREAL, {&cutoff},
473 "Largest distance (nm) to be considered in a cluster" },
474 { "-mol", FALSE, etBOOL, {&bMol},
475 "Cluster molecules rather than atoms (needs [TT].tpr[tt] file)" },
476 { "-pbc", FALSE, etBOOL, {&bPBC},
477 "Use periodic boundary conditions" },
478 { "-nskip", FALSE, etINT, {&nskip},
479 "Number of frames to skip between writing" },
480 { "-nlevels", FALSE, etINT, {&nlevels},
481 "Number of levels of grey in [TT].xpm[tt] output" },
482 { "-ndf", FALSE, etINT, {&ndf},
483 "Number of degrees of freedom of the entire system for temperature calculation. If not set, the number of atoms times three is used." },
484 { "-rgblo", FALSE, etRVEC, {rlo},
485 "RGB values for the color of the lowest occupied cluster size" },
486 { "-rgbhi", FALSE, etRVEC, {rhi},
487 "RGB values for the color of the highest occupied cluster size" }
489 #define NPA asize(pa)
490 const char *fnNDX, *fnTPR;
495 { efTRX, "-f", NULL, ffREAD },
496 { efTPR, NULL, NULL, ffOPTRD },
497 { efNDX, NULL, NULL, ffOPTRD },
498 { efXPM, "-o", "csize", ffWRITE },
499 { efXPM, "-ow", "csizew", ffWRITE },
500 { efXVG, "-nc", "nclust", ffWRITE },
501 { efXVG, "-mc", "maxclust", ffWRITE },
502 { efXVG, "-ac", "avclust", ffWRITE },
503 { efXVG, "-hc", "histo-clust", ffWRITE },
504 { efXVG, "-temp", "temp", ffOPTWR },
505 { efNDX, "-mcn", "maxclust", ffOPTWR }
507 #define NFILE asize(fnm)
509 parse_common_args(&argc, argv,
510 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
511 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
513 fnNDX = ftp2fn_null(efNDX, NFILE, fnm);
514 rgblo.r = rlo[XX], rgblo.g = rlo[YY], rgblo.b = rlo[ZZ];
515 rgbhi.r = rhi[XX], rgbhi.g = rhi[YY], rgbhi.b = rhi[ZZ];
517 fnTPR = ftp2fn_null(efTPR, NFILE, fnm);
520 gmx_fatal(FARGS, "You need a tpr file for the -mol option");
523 clust_size(fnNDX, ftp2fn(efTRX, NFILE, fnm), opt2fn("-o", NFILE, fnm),
524 opt2fn("-ow", NFILE, fnm),
525 opt2fn("-nc", NFILE, fnm), opt2fn("-ac", NFILE, fnm),
526 opt2fn("-mc", NFILE, fnm), opt2fn("-hc", NFILE, fnm),
527 opt2fn("-temp", NFILE, fnm), opt2fn("-mcn", NFILE, fnm),
529 cutoff, nskip, nlevels, rgblo, rgbhi, ndf, oenv);