2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
53 #include "gromacs/utility/smalloc.h"
56 #include "gromacs/math/units.h"
60 #include "gromacs/fileio/matio.h"
61 #include "mtop_util.h"
64 #include "gromacs/utility/fatalerror.h"
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
)
74 FILE *fp
, *gp
, *hp
, *tp
;
75 atom_id
*index
= NULL
;
78 rvec
*x
= NULL
, *v
= NULL
, dx
;
82 gmx_bool bSame
, bTPRwarn
= TRUE
;
86 gmx_mtop_t
*mtop
= NULL
;
89 gmx_mtop_atomlookup_t alook
;
91 int version
, generation
, ii
, jj
, nsame
;
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 };
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)",
109 if (!read_first_frame(oenv
, &status
, trx
, &fr
, TRX_NEED_X
| TRX_READ_V
))
120 read_tpxheader(tpr
, &tpxh
, TRUE
, &version
, &generation
);
121 if (tpxh
.natoms
!= natoms
)
123 gmx_fatal(FARGS
, "tpr (%d atoms) and trajectory (%d atoms) do not match!",
124 tpxh
.natoms
, natoms
);
126 ePBC
= read_tpx(tpr
, NULL
, NULL
, &natoms
, NULL
, NULL
, NULL
, mtop
);
134 tfac
= ndf
/(3.0*natoms
);
141 printf("Using molecules rather than atoms. Not reading index file %s\n",
144 mols
= &(mtop
->mols
);
146 /* Make dummy index */
149 for (i
= 0; (i
< nindex
); i
++)
153 gname
= strdup("mols");
157 rd_index(ndx
, 1, &nindex
, &index
, &gname
);
160 alook
= gmx_mtop_atomlookup_init(mtop
);
162 snew(clust_index
, nindex
);
163 snew(clust_size
, nindex
);
168 for (i
= 0; (i
< nindex
); i
++)
176 if ((nskip
== 0) || ((nskip
> 0) && ((nframe
% nskip
) == 0)))
180 set_pbc(&pbc
, ePBC
, fr
.box
);
185 /* Put all atoms/molecules in their own cluster, with size 1 */
186 for (i
= 0; (i
< nindex
); i
++)
188 /* Cluster index is indexed with atom index number */
190 /* Cluster size is indexed with cluster number */
194 /* Loop over atoms */
195 for (i
= 0; (i
< nindex
); i
++)
200 /* Loop over atoms (only half a matrix) */
201 for (j
= i
+1; (j
< nindex
); j
++)
205 /* If they are not in the same cluster already */
210 /* Compute distance */
214 for (ii
= mols
->index
[ai
]; !bSame
&& (ii
< mols
->index
[ai
+1]); ii
++)
216 for (jj
= mols
->index
[aj
]; !bSame
&& (jj
< mols
->index
[aj
+1]); jj
++)
220 pbc_dx(&pbc
, x
[ii
], x
[jj
], dx
);
224 rvec_sub(x
[ii
], x
[jj
], dx
);
227 bSame
= (dx2
< cut2
);
235 pbc_dx(&pbc
, x
[ai
], x
[aj
], dx
);
239 rvec_sub(x
[ai
], x
[aj
], dx
);
242 bSame
= (dx2
< cut2
);
244 /* If distance less than cut-off */
247 /* Merge clusters: check for all atoms whether they are in
248 * cluster cj and if so, put them in ci
250 for (k
= 0; (k
< nindex
); k
++)
252 if (clust_index
[k
] == cj
)
254 if (clust_size
[cj
] <= 0)
256 gmx_fatal(FARGS
, "negative cluster size %d for element %d",
270 t_x
[n_x
-1] = fr
.time
*tf
;
271 srenew(cs_dist
, n_x
);
272 snew(cs_dist
[n_x
-1], nindex
);
276 for (i
= 0; (i
< nindex
); i
++)
279 if (ci
> max_clust_size
)
287 cs_dist
[n_x
-1][ci
-1] += 1.0;
288 max_size
= max(max_size
, ci
);
296 fprintf(fp
, "%14.6e %10d\n", fr
.time
, nclust
);
299 fprintf(gp
, "%14.6e %10.3f\n", fr
.time
, cav
/nav
);
301 fprintf(hp
, "%14.6e %10d\n", fr
.time
, max_clust_size
);
303 /* Analyse velocities, if present */
310 printf("You need a [TT].tpr[tt] file to analyse temperatures\n");
317 /* Loop over clusters and for each cluster compute 1/2 m v^2 */
318 if (max_clust_ind
>= 0)
321 for (i
= 0; (i
< nindex
); i
++)
323 if (clust_index
[i
] == max_clust_ind
)
326 gmx_mtop_atomnr_to_atom(alook
, ai
, &atom
);
327 ekin
+= 0.5*atom
->m
*iprod(v
[ai
], v
[ai
]);
330 temp
= (ekin
*2.0)/(3.0*tfac
*max_clust_size
*BOLTZ
);
331 fprintf(tp
, "%10.3f %10.3f\n", fr
.time
, temp
);
337 while (read_next_frame(oenv
, status
, &fr
));
344 gmx_mtop_atomlookup_destroy(alook
);
346 if (max_clust_ind
>= 0)
348 fp
= gmx_ffopen(mcn
, "w");
349 fprintf(fp
, "[ max_clust ]\n");
350 for (i
= 0; (i
< nindex
); i
++)
352 if (clust_index
[i
] == max_clust_ind
)
356 for (j
= mols
->index
[i
]; (j
< mols
->index
[i
+1]); j
++)
358 fprintf(fp
, "%d\n", j
+1);
363 fprintf(fp
, "%d\n", index
[i
]+1);
370 /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
371 fp
= xvgropen(histo
, "Cluster size distribution", "Cluster size", "()", oenv
);
373 fprintf(fp
, "%5d %8.3f\n", 0, 0.0);
374 for (j
= 0; (j
< max_size
); j
++)
377 for (i
= 0; (i
< n_x
); i
++)
379 nelem
+= cs_dist
[i
][j
];
381 fprintf(fp
, "%5d %8.3f\n", j
+1, nelem
/n_x
);
382 nhisto
+= (int)((j
+1)*nelem
/n_x
);
384 fprintf(fp
, "%5d %8.3f\n", j
+1, 0.0);
387 fprintf(stderr
, "Total number of atoms in clusters = %d\n", nhisto
);
389 /* Look for the smallest entry that is not zero
390 * This will make that zero is white, and not zero is coloured.
394 for (i
= 0; (i
< n_x
); i
++)
396 for (j
= 0; (j
< max_size
); j
++)
398 if ((cs_dist
[i
][j
] > 0) && (cs_dist
[i
][j
] < cmid
))
400 cmid
= cs_dist
[i
][j
];
402 cmax
= max(cs_dist
[i
][j
], cmax
);
405 fprintf(stderr
, "cmid: %g, cmax: %g, max_size: %d\n", cmid
, cmax
, max_size
);
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
);
414 for (i
= 0; (i
< n_x
); i
++)
416 for (j
= 0; (j
< max_size
); j
++)
418 cs_dist
[i
][j
] *= (j
+1);
419 if ((cs_dist
[i
][j
] > 0) && (cs_dist
[i
][j
] < cmid
))
421 cmid
= cs_dist
[i
][j
];
423 cmax
= max(cs_dist
[i
][j
], cmax
);
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
);
438 int gmx_clustsize(int argc
, char *argv
[])
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."
459 static real cutoff
= 0.35;
460 static int nskip
= 0;
461 static int nlevels
= 20;
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 };
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" }
488 #define NPA asize(pa)
489 const char *fnNDX
, *fnTPR
;
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
}
506 #define NFILE asize(fnm)
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
))
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
];
519 fnTPR
= ftp2fn_null(efTPR
, NFILE
, fnm
);
522 gmx_fatal(FARGS
, "You need a tpr file for the -mol option");
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
),
531 cutoff
, nskip
, nlevels
, rgblo
, rgbhi
, ndf
, oenv
);