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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
67 #include "mtop_util.h"
71 static void clust_size(const char *ndx,const char *trx,const char *xpm,
72 const char *xpmw,const char *ncl,const char *acl,
73 const char *mcl,const char *histo,const char *tempf,
74 const char *mcn,gmx_bool bMol,gmx_bool bPBC,const char *tpr,
75 real cut,int nskip,int nlevels,
76 t_rgb rmid,t_rgb rhi,int ndf,
77 const output_env_t oenv)
83 rvec *x=NULL,*v=NULL,dx;
87 gmx_bool bSame,bTPRwarn=TRUE;
91 gmx_mtop_t *mtop=NULL;
94 gmx_mtop_atomlookup_t alook;
96 int version,generation,ii,jj,nsame;
98 /* Cluster size distribution (matrix) */
100 real tf,dx2,cut2,*t_x=NULL,*t_y,cmid,cmax,cav,ekin;
101 int i,j,k,ai,aj,ak,ci,cj,nframe,nclust,n_x,n_y,max_size=0;
102 int *clust_index,*clust_size,max_clust_size,max_clust_ind,nav,nhisto;
103 t_rgb rlo = { 1.0, 1.0, 1.0 };
105 clear_trxframe(&fr,TRUE);
106 sprintf(timebuf,"Time (%s)",output_env_get_time_unit(oenv));
107 tf = output_env_get_time_factor(oenv);
108 fp = xvgropen(ncl,"Number of clusters",timebuf,"N",oenv);
109 gp = xvgropen(acl,"Average cluster size",timebuf,"#molecules",oenv);
110 hp = xvgropen(mcl,"Max cluster size",timebuf,"#molecules",oenv);
111 tp = xvgropen(tempf,"Temperature of largest cluster",timebuf,"T (K)",
114 if (!read_first_frame(oenv,&status,trx,&fr,TRX_NEED_X | TRX_READ_V))
122 read_tpxheader(tpr,&tpxh,TRUE,&version,&generation);
123 if (tpxh.natoms != natoms)
124 gmx_fatal(FARGS,"tpr (%d atoms) and xtc (%d atoms) do not match!",
126 ePBC = read_tpx(tpr,NULL,NULL,&natoms,NULL,NULL,NULL,mtop);
131 tfac = ndf/(3.0*natoms);
135 printf("Using molecules rather than atoms. Not reading index file %s\n",
137 mols = &(mtop->mols);
139 /* Make dummy index */
142 for(i=0; (i<nindex); i++)
144 gname = strdup("mols");
147 rd_index(ndx,1,&nindex,&index,&gname);
149 alook = gmx_mtop_atomlookup_init(mtop);
151 snew(clust_index,nindex);
152 snew(clust_size,nindex);
157 for(i=0; (i<nindex); i++)
162 if ((nskip == 0) || ((nskip > 0) && ((nframe % nskip) == 0))) {
164 set_pbc(&pbc,ePBC,fr.box);
168 /* Put all atoms/molecules in their own cluster, with size 1 */
169 for(i=0; (i<nindex); i++) {
170 /* Cluster index is indexed with atom index number */
172 /* Cluster size is indexed with cluster number */
176 /* Loop over atoms */
177 for(i=0; (i<nindex); i++) {
181 /* Loop over atoms (only half a matrix) */
182 for(j=i+1; (j<nindex); j++) {
185 /* If they are not in the same cluster already */
189 /* Compute distance */
192 for(ii=mols->index[ai]; !bSame && (ii<mols->index[ai+1]); ii++) {
193 for(jj=mols->index[aj]; !bSame && (jj<mols->index[aj+1]); jj++) {
195 pbc_dx(&pbc,x[ii],x[jj],dx);
197 rvec_sub(x[ii],x[jj],dx);
199 bSame = (dx2 < cut2);
205 pbc_dx(&pbc,x[ai],x[aj],dx);
207 rvec_sub(x[ai],x[aj],dx);
209 bSame = (dx2 < cut2);
211 /* If distance less than cut-off */
213 /* Merge clusters: check for all atoms whether they are in
214 * cluster cj and if so, put them in ci
216 for(k=0; (k<nindex); k++) {
217 if (clust_index[k] == cj) {
218 if (clust_size[cj] <= 0)
219 gmx_fatal(FARGS,"negative cluster size %d for element %d",
232 t_x[n_x-1] = fr.time*tf;
234 snew(cs_dist[n_x-1],nindex);
238 for(i=0; (i<nindex); i++) {
240 if (ci > max_clust_size) {
246 cs_dist[n_x-1][ci-1] += 1.0;
247 max_size = max(max_size,ci);
254 fprintf(fp,"%14.6e %10d\n",fr.time,nclust);
256 fprintf(gp,"%14.6e %10.3f\n",fr.time,cav/nav);
257 fprintf(hp, "%14.6e %10d\n",fr.time,max_clust_size);
259 /* Analyse velocities, if present */
263 printf("You need a [TT].tpr[tt] file to analyse temperatures\n");
269 /* Loop over clusters and for each cluster compute 1/2 m v^2 */
270 if (max_clust_ind >= 0) {
272 for(i=0; (i<nindex); i++)
273 if (clust_index[i] == max_clust_ind) {
275 gmx_mtop_atomnr_to_atom(alook,ai,&atom);
276 ekin += 0.5*atom->m*iprod(v[ai],v[ai]);
278 temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
279 fprintf(tp,"%10.3f %10.3f\n",fr.time,temp);
284 } while (read_next_frame(oenv,status,&fr));
291 gmx_mtop_atomlookup_destroy(alook);
293 if (max_clust_ind >= 0) {
294 fp = ffopen(mcn,"w");
295 fprintf(fp,"[ max_clust ]\n");
296 for(i=0; (i<nindex); i++)
297 if (clust_index[i] == max_clust_ind) {
299 for(j=mols->index[i]; (j<mols->index[i+1]); j++)
300 fprintf(fp,"%d\n",j+1);
303 fprintf(fp,"%d\n",index[i]+1);
309 /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
310 fp = xvgropen(histo,"Cluster size distribution","Cluster size","()",oenv);
312 fprintf(fp,"%5d %8.3f\n",0,0.0);
313 for(j=0; (j<max_size); j++) {
315 for(i=0; (i<n_x); i++)
316 nelem += cs_dist[i][j];
317 fprintf(fp,"%5d %8.3f\n",j+1,nelem/n_x);
318 nhisto += (int)((j+1)*nelem/n_x);
320 fprintf(fp,"%5d %8.3f\n",j+1,0.0);
323 fprintf(stderr,"Total number of atoms in clusters = %d\n",nhisto);
325 /* Look for the smallest entry that is not zero
326 * This will make that zero is white, and not zero is coloured.
330 for(i=0; (i<n_x); i++)
331 for(j=0; (j<max_size); j++) {
332 if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
333 cmid = cs_dist[i][j];
334 cmax = max(cs_dist[i][j],cmax);
336 fprintf(stderr,"cmid: %g, cmax: %g, max_size: %d\n",cmid,cmax,max_size);
338 fp = ffopen(xpm,"w");
339 write_xpm3(fp,0,"Cluster size distribution","# clusters",timebuf,"Size",
340 n_x,max_size,t_x,t_y,cs_dist,0,cmid,cmax,
341 rlo,rmid,rhi,&nlevels);
345 for(i=0; (i<n_x); i++)
346 for(j=0; (j<max_size); j++) {
347 cs_dist[i][j] *= (j+1);
348 if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
349 cmid = cs_dist[i][j];
350 cmax = max(cs_dist[i][j],cmax);
352 fprintf(stderr,"cmid: %g, cmax: %g, max_size: %d\n",cmid,cmax,max_size);
353 fp = ffopen(xpmw,"w");
354 write_xpm3(fp,0,"Weighted cluster size distribution","Fraction",timebuf,
355 "Size", n_x,max_size,t_x,t_y,cs_dist,0,cmid,cmax,
356 rlo,rmid,rhi,&nlevels);
364 int gmx_clustsize(int argc,char *argv[])
366 const char *desc[] = {
367 "This program computes the size distributions of molecular/atomic clusters in",
368 "the gas phase. The output is given in the form of an [TT].xpm[tt] file.",
369 "The total number of clusters is written to an [TT].xvg[tt] file.[PAR]",
370 "When the [TT]-mol[tt] option is given clusters will be made out of",
371 "molecules rather than atoms, which allows clustering of large molecules.",
372 "In this case an index file would still contain atom numbers",
373 "or your calculation will die with a SEGV.[PAR]",
374 "When velocities are present in your trajectory, the temperature of",
375 "the largest cluster will be printed in a separate [TT].xvg[tt] file assuming",
376 "that the particles are free to move. If you are using constraints,",
377 "please correct the temperature. For instance water simulated with SHAKE",
378 "or SETTLE will yield a temperature that is 1.5 times too low. You can",
379 "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal",
380 "of center of mass motion into account.[PAR]",
381 "The [TT]-mc[tt] option will produce an index file containing the",
382 "atom numbers of the largest cluster."
385 static real cutoff = 0.35;
386 static int nskip = 0;
387 static int nlevels = 20;
389 static gmx_bool bMol = FALSE;
390 static gmx_bool bPBC = TRUE;
391 static rvec rlo = { 1.0, 1.0, 0.0 };
392 static rvec rhi = { 0.0, 0.0, 1.0 };
397 { "-cut", FALSE, etREAL, {&cutoff},
398 "Largest distance (nm) to be considered in a cluster" },
399 { "-mol", FALSE, etBOOL, {&bMol},
400 "Cluster molecules rather than atoms (needs [TT].tpr[tt] file)" },
401 { "-pbc", FALSE, etBOOL, {&bPBC},
402 "Use periodic boundary conditions" },
403 { "-nskip", FALSE, etINT, {&nskip},
404 "Number of frames to skip between writing" },
405 { "-nlevels", FALSE, etINT, {&nlevels},
406 "Number of levels of grey in [TT].xpm[tt] output" },
407 { "-ndf", FALSE, etINT, {&ndf},
408 "Number of degrees of freedom of the entire system for temperature calculation. If not set, the number of atoms times three is used." },
409 { "-rgblo", FALSE, etRVEC, {rlo},
410 "RGB values for the color of the lowest occupied cluster size" },
411 { "-rgbhi", FALSE, etRVEC, {rhi},
412 "RGB values for the color of the highest occupied cluster size" }
414 #define NPA asize(pa)
415 const char *fnNDX,*fnTPR;
420 { efTRX, "-f", NULL, ffREAD },
421 { efTPR, NULL, NULL, ffOPTRD },
422 { efNDX, NULL, NULL, ffOPTRD },
423 { efXPM, "-o", "csize", ffWRITE },
424 { efXPM, "-ow","csizew", ffWRITE },
425 { efXVG, "-nc","nclust", ffWRITE },
426 { efXVG, "-mc","maxclust", ffWRITE },
427 { efXVG, "-ac","avclust", ffWRITE },
428 { efXVG, "-hc","histo-clust", ffWRITE },
429 { efXVG, "-temp","temp", ffOPTWR },
430 { efNDX, "-mcn", "maxclust", ffOPTWR }
432 #define NFILE asize(fnm)
434 CopyRight(stderr,argv[0]);
435 parse_common_args(&argc,argv,
436 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
437 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
439 fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
440 rgblo.r = rlo[XX],rgblo.g = rlo[YY],rgblo.b = rlo[ZZ];
441 rgbhi.r = rhi[XX],rgbhi.g = rhi[YY],rgbhi.b = rhi[ZZ];
443 fnTPR = ftp2fn_null(efTPR,NFILE,fnm);
445 gmx_fatal(FARGS,"You need a tpr file for the -mol option");
447 clust_size(fnNDX,ftp2fn(efTRX,NFILE,fnm),opt2fn("-o",NFILE,fnm),
448 opt2fn("-ow",NFILE,fnm),
449 opt2fn("-nc",NFILE,fnm),opt2fn("-ac",NFILE,fnm),
450 opt2fn("-mc",NFILE,fnm),opt2fn("-hc",NFILE,fnm),
451 opt2fn("-temp",NFILE,fnm),opt2fn("-mcn",NFILE,fnm),
453 cutoff,nskip,nlevels,rgblo,rgbhi,ndf,oenv);