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
64 #include "mtop_util.h"
68 static void clust_size(const char *ndx
,const char *trx
,const char *xpm
,
69 const char *xpmw
,const char *ncl
,const char *acl
,
70 const char *mcl
,const char *histo
,const char *tempf
,
71 const char *mcn
,gmx_bool bMol
,gmx_bool bPBC
,const char *tpr
,
72 real cut
,int nskip
,int nlevels
,
73 t_rgb rmid
,t_rgb rhi
,int ndf
,
74 const output_env_t oenv
)
80 rvec
*x
=NULL
,*v
=NULL
,dx
;
84 gmx_bool bSame
,bTPRwarn
=TRUE
;
88 gmx_mtop_t
*mtop
=NULL
;
91 gmx_mtop_atomlookup_t alook
;
93 int version
,generation
,ii
,jj
,nsame
;
95 /* Cluster size distribution (matrix) */
97 real tf
,dx2
,cut2
,*t_x
=NULL
,*t_y
,cmid
,cmax
,cav
,ekin
;
98 int i
,j
,k
,ai
,aj
,ak
,ci
,cj
,nframe
,nclust
,n_x
,n_y
,max_size
=0;
99 int *clust_index
,*clust_size
,max_clust_size
,max_clust_ind
,nav
,nhisto
;
100 t_rgb rlo
= { 1.0, 1.0, 1.0 };
102 clear_trxframe(&fr
,TRUE
);
103 sprintf(timebuf
,"Time (%s)",output_env_get_time_unit(oenv
));
104 tf
= output_env_get_time_factor(oenv
);
105 fp
= xvgropen(ncl
,"Number of clusters",timebuf
,"N",oenv
);
106 gp
= xvgropen(acl
,"Average cluster size",timebuf
,"#molecules",oenv
);
107 hp
= xvgropen(mcl
,"Max cluster size",timebuf
,"#molecules",oenv
);
108 tp
= xvgropen(tempf
,"Temperature of largest cluster",timebuf
,"T (K)",
111 if (!read_first_frame(oenv
,&status
,trx
,&fr
,TRX_NEED_X
| TRX_READ_V
))
119 read_tpxheader(tpr
,&tpxh
,TRUE
,&version
,&generation
);
120 if (tpxh
.natoms
!= natoms
)
121 gmx_fatal(FARGS
,"tpr (%d atoms) and xtc (%d atoms) do not match!",
123 ePBC
= read_tpx(tpr
,NULL
,NULL
,&natoms
,NULL
,NULL
,NULL
,mtop
);
128 tfac
= ndf
/(3.0*natoms
);
132 printf("Using molecules rather than atoms. Not reading index file %s\n",
134 mols
= &(mtop
->mols
);
136 /* Make dummy index */
139 for(i
=0; (i
<nindex
); i
++)
141 gname
= strdup("mols");
144 rd_index(ndx
,1,&nindex
,&index
,&gname
);
146 alook
= gmx_mtop_atomlookup_init(mtop
);
148 snew(clust_index
,nindex
);
149 snew(clust_size
,nindex
);
154 for(i
=0; (i
<nindex
); i
++)
159 if ((nskip
== 0) || ((nskip
> 0) && ((nframe
% nskip
) == 0))) {
161 set_pbc(&pbc
,ePBC
,fr
.box
);
165 /* Put all atoms/molecules in their own cluster, with size 1 */
166 for(i
=0; (i
<nindex
); i
++) {
167 /* Cluster index is indexed with atom index number */
169 /* Cluster size is indexed with cluster number */
173 /* Loop over atoms */
174 for(i
=0; (i
<nindex
); i
++) {
178 /* Loop over atoms (only half a matrix) */
179 for(j
=i
+1; (j
<nindex
); j
++) {
182 /* If they are not in the same cluster already */
186 /* Compute distance */
189 for(ii
=mols
->index
[ai
]; !bSame
&& (ii
<mols
->index
[ai
+1]); ii
++) {
190 for(jj
=mols
->index
[aj
]; !bSame
&& (jj
<mols
->index
[aj
+1]); jj
++) {
192 pbc_dx(&pbc
,x
[ii
],x
[jj
],dx
);
194 rvec_sub(x
[ii
],x
[jj
],dx
);
196 bSame
= (dx2
< cut2
);
202 pbc_dx(&pbc
,x
[ai
],x
[aj
],dx
);
204 rvec_sub(x
[ai
],x
[aj
],dx
);
206 bSame
= (dx2
< cut2
);
208 /* If distance less than cut-off */
210 /* Merge clusters: check for all atoms whether they are in
211 * cluster cj and if so, put them in ci
213 for(k
=0; (k
<nindex
); k
++) {
214 if (clust_index
[k
] == cj
) {
215 if (clust_size
[cj
] <= 0)
216 gmx_fatal(FARGS
,"negative cluster size %d for element %d",
229 t_x
[n_x
-1] = fr
.time
*tf
;
231 snew(cs_dist
[n_x
-1],nindex
);
235 for(i
=0; (i
<nindex
); i
++) {
237 if (ci
> max_clust_size
) {
243 cs_dist
[n_x
-1][ci
-1] += 1.0;
244 max_size
= max(max_size
,ci
);
251 fprintf(fp
,"%14.6e %10d\n",fr
.time
,nclust
);
253 fprintf(gp
,"%14.6e %10.3f\n",fr
.time
,cav
/nav
);
254 fprintf(hp
, "%14.6e %10d\n",fr
.time
,max_clust_size
);
256 /* Analyse velocities, if present */
260 printf("You need a [TT].tpr[tt] file to analyse temperatures\n");
266 /* Loop over clusters and for each cluster compute 1/2 m v^2 */
267 if (max_clust_ind
>= 0) {
269 for(i
=0; (i
<nindex
); i
++)
270 if (clust_index
[i
] == max_clust_ind
) {
272 gmx_mtop_atomnr_to_atom(alook
,ai
,&atom
);
273 ekin
+= 0.5*atom
->m
*iprod(v
[ai
],v
[ai
]);
275 temp
= (ekin
*2.0)/(3.0*tfac
*max_clust_size
*BOLTZ
);
276 fprintf(tp
,"%10.3f %10.3f\n",fr
.time
,temp
);
281 } while (read_next_frame(oenv
,status
,&fr
));
288 gmx_mtop_atomlookup_destroy(alook
);
290 if (max_clust_ind
>= 0) {
291 fp
= ffopen(mcn
,"w");
292 fprintf(fp
,"[ max_clust ]\n");
293 for(i
=0; (i
<nindex
); i
++)
294 if (clust_index
[i
] == max_clust_ind
) {
296 for(j
=mols
->index
[i
]; (j
<mols
->index
[i
+1]); j
++)
297 fprintf(fp
,"%d\n",j
+1);
300 fprintf(fp
,"%d\n",index
[i
]+1);
306 /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
307 fp
= xvgropen(histo
,"Cluster size distribution","Cluster size","()",oenv
);
309 fprintf(fp
,"%5d %8.3f\n",0,0.0);
310 for(j
=0; (j
<max_size
); j
++) {
312 for(i
=0; (i
<n_x
); i
++)
313 nelem
+= cs_dist
[i
][j
];
314 fprintf(fp
,"%5d %8.3f\n",j
+1,nelem
/n_x
);
315 nhisto
+= (int)((j
+1)*nelem
/n_x
);
317 fprintf(fp
,"%5d %8.3f\n",j
+1,0.0);
320 fprintf(stderr
,"Total number of atoms in clusters = %d\n",nhisto
);
322 /* Look for the smallest entry that is not zero
323 * This will make that zero is white, and not zero is coloured.
327 for(i
=0; (i
<n_x
); i
++)
328 for(j
=0; (j
<max_size
); j
++) {
329 if ((cs_dist
[i
][j
] > 0) && (cs_dist
[i
][j
] < cmid
))
330 cmid
= cs_dist
[i
][j
];
331 cmax
= max(cs_dist
[i
][j
],cmax
);
333 fprintf(stderr
,"cmid: %g, cmax: %g, max_size: %d\n",cmid
,cmax
,max_size
);
335 fp
= ffopen(xpm
,"w");
336 write_xpm3(fp
,0,"Cluster size distribution","# clusters",timebuf
,"Size",
337 n_x
,max_size
,t_x
,t_y
,cs_dist
,0,cmid
,cmax
,
338 rlo
,rmid
,rhi
,&nlevels
);
342 for(i
=0; (i
<n_x
); i
++)
343 for(j
=0; (j
<max_size
); j
++) {
344 cs_dist
[i
][j
] *= (j
+1);
345 if ((cs_dist
[i
][j
] > 0) && (cs_dist
[i
][j
] < cmid
))
346 cmid
= cs_dist
[i
][j
];
347 cmax
= max(cs_dist
[i
][j
],cmax
);
349 fprintf(stderr
,"cmid: %g, cmax: %g, max_size: %d\n",cmid
,cmax
,max_size
);
350 fp
= ffopen(xpmw
,"w");
351 write_xpm3(fp
,0,"Weighted cluster size distribution","Fraction",timebuf
,
352 "Size", n_x
,max_size
,t_x
,t_y
,cs_dist
,0,cmid
,cmax
,
353 rlo
,rmid
,rhi
,&nlevels
);
361 int gmx_clustsize(int argc
,char *argv
[])
363 const char *desc
[] = {
364 "This program computes the size distributions of molecular/atomic clusters in",
365 "the gas phase. The output is given in the form of an [TT].xpm[tt] file.",
366 "The total number of clusters is written to an [TT].xvg[tt] file.[PAR]",
367 "When the [TT]-mol[tt] option is given clusters will be made out of",
368 "molecules rather than atoms, which allows clustering of large molecules.",
369 "In this case an index file would still contain atom numbers",
370 "or your calculation will die with a SEGV.[PAR]",
371 "When velocities are present in your trajectory, the temperature of",
372 "the largest cluster will be printed in a separate [TT].xvg[tt] file assuming",
373 "that the particles are free to move. If you are using constraints,",
374 "please correct the temperature. For instance water simulated with SHAKE",
375 "or SETTLE will yield a temperature that is 1.5 times too low. You can",
376 "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal",
377 "of center of mass motion into account.[PAR]",
378 "The [TT]-mc[tt] option will produce an index file containing the",
379 "atom numbers of the largest cluster."
382 static real cutoff
= 0.35;
383 static int nskip
= 0;
384 static int nlevels
= 20;
386 static gmx_bool bMol
= FALSE
;
387 static gmx_bool bPBC
= TRUE
;
388 static rvec rlo
= { 1.0, 1.0, 0.0 };
389 static rvec rhi
= { 0.0, 0.0, 1.0 };
394 { "-cut", FALSE
, etREAL
, {&cutoff
},
395 "Largest distance (nm) to be considered in a cluster" },
396 { "-mol", FALSE
, etBOOL
, {&bMol
},
397 "Cluster molecules rather than atoms (needs [TT].tpr[tt] file)" },
398 { "-pbc", FALSE
, etBOOL
, {&bPBC
},
399 "Use periodic boundary conditions" },
400 { "-nskip", FALSE
, etINT
, {&nskip
},
401 "Number of frames to skip between writing" },
402 { "-nlevels", FALSE
, etINT
, {&nlevels
},
403 "Number of levels of grey in [TT].xpm[tt] output" },
404 { "-ndf", FALSE
, etINT
, {&ndf
},
405 "Number of degrees of freedom of the entire system for temperature calculation. If not set, the number of atoms times three is used." },
406 { "-rgblo", FALSE
, etRVEC
, {rlo
},
407 "RGB values for the color of the lowest occupied cluster size" },
408 { "-rgbhi", FALSE
, etRVEC
, {rhi
},
409 "RGB values for the color of the highest occupied cluster size" }
411 #define NPA asize(pa)
412 const char *fnNDX
,*fnTPR
;
417 { efTRX
, "-f", NULL
, ffREAD
},
418 { efTPR
, NULL
, NULL
, ffOPTRD
},
419 { efNDX
, NULL
, NULL
, ffOPTRD
},
420 { efXPM
, "-o", "csize", ffWRITE
},
421 { efXPM
, "-ow","csizew", ffWRITE
},
422 { efXVG
, "-nc","nclust", ffWRITE
},
423 { efXVG
, "-mc","maxclust", ffWRITE
},
424 { efXVG
, "-ac","avclust", ffWRITE
},
425 { efXVG
, "-hc","histo-clust", ffWRITE
},
426 { efXVG
, "-temp","temp", ffOPTWR
},
427 { efNDX
, "-mcn", "maxclust", ffOPTWR
}
429 #define NFILE asize(fnm)
431 CopyRight(stderr
,argv
[0]);
432 parse_common_args(&argc
,argv
,
433 PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_TIME_UNIT
| PCA_BE_NICE
,
434 NFILE
,fnm
,NPA
,pa
,asize(desc
),desc
,0,NULL
,&oenv
);
436 fnNDX
= ftp2fn_null(efNDX
,NFILE
,fnm
);
437 rgblo
.r
= rlo
[XX
],rgblo
.g
= rlo
[YY
],rgblo
.b
= rlo
[ZZ
];
438 rgbhi
.r
= rhi
[XX
],rgbhi
.g
= rhi
[YY
],rgbhi
.b
= rhi
[ZZ
];
440 fnTPR
= ftp2fn_null(efTPR
,NFILE
,fnm
);
442 gmx_fatal(FARGS
,"You need a tpr file for the -mol option");
444 clust_size(fnNDX
,ftp2fn(efTRX
,NFILE
,fnm
),opt2fn("-o",NFILE
,fnm
),
445 opt2fn("-ow",NFILE
,fnm
),
446 opt2fn("-nc",NFILE
,fnm
),opt2fn("-ac",NFILE
,fnm
),
447 opt2fn("-mc",NFILE
,fnm
),opt2fn("-hc",NFILE
,fnm
),
448 opt2fn("-temp",NFILE
,fnm
),opt2fn("-mcn",NFILE
,fnm
),
450 cutoff
,nskip
,nlevels
,rgblo
,rgbhi
,ndf
,oenv
);