/*
- *
+ *
* This source code is part of
- *
+ *
* G R O M A C S
- *
+ *
* GROningen MAchine for Chemical Simulations
- *
+ *
* VERSION 3.3.2
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
- *
+ *
* If you want to redistribute modifications, please consider that
* scientific software is very special. Version control is crucial -
* bugs must be traceable. We will be happy to consider code for
* inclusion in the official distribution, but derived work must not
* be called official GROMACS. Details are found in the README & COPYING
* files - if they are missing, get the official version at www.gromacs.org.
- *
+ *
* To help us fund GROMACS development, we humbly ask that you cite
* the papers on the package - you can find them in the top README file.
- *
+ *
* For more info, check our website at http://www.gromacs.org
- *
+ *
* And Hey:
* Groningen Machine for Chemical Simulation
*/
#include "gmx_ana.h"
-static void clust_size(const char *ndx,const char *trx,const char *xpm,
- const char *xpmw,const char *ncl,const char *acl,
- const char *mcl,const char *histo,const char *tempf,
- const char *mcn,gmx_bool bMol,gmx_bool bPBC,const char *tpr,
- real cut,int nskip,int nlevels,
- t_rgb rmid,t_rgb rhi,int ndf,
+static void clust_size(const char *ndx, const char *trx, const char *xpm,
+ const char *xpmw, const char *ncl, const char *acl,
+ const char *mcl, const char *histo, const char *tempf,
+ const char *mcn, gmx_bool bMol, gmx_bool bPBC, const char *tpr,
+ real cut, int nskip, int nlevels,
+ t_rgb rmid, t_rgb rhi, int ndf,
const output_env_t oenv)
{
- FILE *fp,*gp,*hp,*tp;
- atom_id *index=NULL;
- int nindex,natoms;
- t_trxstatus *status;
- rvec *x=NULL,*v=NULL,dx;
- t_pbc pbc;
- char *gname;
- char timebuf[32];
- gmx_bool bSame,bTPRwarn=TRUE;
- /* Topology stuff */
- t_trxframe fr;
- t_tpxheader tpxh;
- gmx_mtop_t *mtop=NULL;
- int ePBC=-1;
- t_block *mols=NULL;
- gmx_mtop_atomlookup_t alook;
- t_atom *atom;
- int version,generation,ii,jj,nsame;
- real temp,tfac;
- /* Cluster size distribution (matrix) */
- real **cs_dist=NULL;
- real tf,dx2,cut2,*t_x=NULL,*t_y,cmid,cmax,cav,ekin;
- int i,j,k,ai,aj,ak,ci,cj,nframe,nclust,n_x,n_y,max_size=0;
- int *clust_index,*clust_size,max_clust_size,max_clust_ind,nav,nhisto;
- t_rgb rlo = { 1.0, 1.0, 1.0 };
-
- clear_trxframe(&fr,TRUE);
- sprintf(timebuf,"Time (%s)",output_env_get_time_unit(oenv));
- tf = output_env_get_time_factor(oenv);
- fp = xvgropen(ncl,"Number of clusters",timebuf,"N",oenv);
- gp = xvgropen(acl,"Average cluster size",timebuf,"#molecules",oenv);
- hp = xvgropen(mcl,"Max cluster size",timebuf,"#molecules",oenv);
- tp = xvgropen(tempf,"Temperature of largest cluster",timebuf,"T (K)",
- oenv);
-
- if (!read_first_frame(oenv,&status,trx,&fr,TRX_NEED_X | TRX_READ_V))
- gmx_file(trx);
-
- natoms = fr.natoms;
- x = fr.x;
-
- if (tpr) {
- snew(mtop,1);
- read_tpxheader(tpr,&tpxh,TRUE,&version,&generation);
- if (tpxh.natoms != natoms)
- gmx_fatal(FARGS,"tpr (%d atoms) and xtc (%d atoms) do not match!",
- tpxh.natoms,natoms);
- ePBC = read_tpx(tpr,NULL,NULL,&natoms,NULL,NULL,NULL,mtop);
- }
- if (ndf <= -1)
- tfac = 1;
- else
- tfac = ndf/(3.0*natoms);
-
- if (bMol) {
- if (ndx)
- printf("Using molecules rather than atoms. Not reading index file %s\n",
- ndx);
- mols = &(mtop->mols);
-
- /* Make dummy index */
- nindex = mols->nr;
- snew(index,nindex);
- for(i=0; (i<nindex); i++)
- index[i] = i;
- gname = strdup("mols");
- }
- else
- rd_index(ndx,1,&nindex,&index,&gname);
-
- alook = gmx_mtop_atomlookup_init(mtop);
-
- snew(clust_index,nindex);
- snew(clust_size,nindex);
- cut2 = cut*cut;
- nframe = 0;
- n_x = 0;
- snew(t_y,nindex);
- for(i=0; (i<nindex); i++)
- t_y[i] = i+1;
- max_clust_size = 1;
- max_clust_ind = -1;
- do {
- if ((nskip == 0) || ((nskip > 0) && ((nframe % nskip) == 0))) {
- if (bPBC)
- set_pbc(&pbc,ePBC,fr.box);
- max_clust_size = 1;
- max_clust_ind = -1;
-
- /* Put all atoms/molecules in their own cluster, with size 1 */
- for(i=0; (i<nindex); i++) {
- /* Cluster index is indexed with atom index number */
- clust_index[i] = i;
- /* Cluster size is indexed with cluster number */
- clust_size[i] = 1;
- }
-
- /* Loop over atoms */
- for(i=0; (i<nindex); i++) {
- ai = index[i];
- ci = clust_index[i];
-
- /* Loop over atoms (only half a matrix) */
- for(j=i+1; (j<nindex); j++) {
- cj = clust_index[j];
-
- /* If they are not in the same cluster already */
- if (ci != cj) {
- aj = index[j];
-
- /* Compute distance */
- if (bMol) {
- bSame = FALSE;
- for(ii=mols->index[ai]; !bSame && (ii<mols->index[ai+1]); ii++) {
- for(jj=mols->index[aj]; !bSame && (jj<mols->index[aj+1]); jj++) {
- if (bPBC)
- pbc_dx(&pbc,x[ii],x[jj],dx);
- else
- rvec_sub(x[ii],x[jj],dx);
- dx2 = iprod(dx,dx);
- bSame = (dx2 < cut2);
- }
- }
- }
- else {
- if (bPBC)
- pbc_dx(&pbc,x[ai],x[aj],dx);
- else
- rvec_sub(x[ai],x[aj],dx);
- dx2 = iprod(dx,dx);
- bSame = (dx2 < cut2);
- }
- /* If distance less than cut-off */
- if (bSame) {
- /* Merge clusters: check for all atoms whether they are in
- * cluster cj and if so, put them in ci
- */
- for(k=0; (k<nindex); k++) {
- if (clust_index[k] == cj) {
- if (clust_size[cj] <= 0)
- gmx_fatal(FARGS,"negative cluster size %d for element %d",
- clust_size[cj],cj);
- clust_size[cj]--;
- clust_index[k] = ci;
- clust_size[ci]++;
- }
- }
- }
- }
- }
- }
- n_x++;
- srenew(t_x,n_x);
- t_x[n_x-1] = fr.time*tf;
- srenew(cs_dist,n_x);
- snew(cs_dist[n_x-1],nindex);
- nclust = 0;
- cav = 0;
- nav = 0;
- for(i=0; (i<nindex); i++) {
- ci = clust_size[i];
- if (ci > max_clust_size) {
- max_clust_size = ci;
- max_clust_ind = i;
- }
- if (ci > 0) {
- nclust++;
- cs_dist[n_x-1][ci-1] += 1.0;
- max_size = max(max_size,ci);
- if (ci > 1) {
- cav += ci;
- nav++;
- }
- }
- }
- fprintf(fp,"%14.6e %10d\n",fr.time,nclust);
- if (nav > 0)
- fprintf(gp,"%14.6e %10.3f\n",fr.time,cav/nav);
- fprintf(hp, "%14.6e %10d\n",fr.time,max_clust_size);
+ FILE *fp, *gp, *hp, *tp;
+ atom_id *index = NULL;
+ int nindex, natoms;
+ t_trxstatus *status;
+ rvec *x = NULL, *v = NULL, dx;
+ t_pbc pbc;
+ char *gname;
+ char timebuf[32];
+ gmx_bool bSame, bTPRwarn = TRUE;
+ /* Topology stuff */
+ t_trxframe fr;
+ t_tpxheader tpxh;
+ gmx_mtop_t *mtop = NULL;
+ int ePBC = -1;
+ t_block *mols = NULL;
+ gmx_mtop_atomlookup_t alook;
+ t_atom *atom;
+ int version, generation, ii, jj, nsame;
+ real temp, tfac;
+ /* Cluster size distribution (matrix) */
+ real **cs_dist = NULL;
+ real tf, dx2, cut2, *t_x = NULL, *t_y, cmid, cmax, cav, ekin;
+ int i, j, k, ai, aj, ak, ci, cj, nframe, nclust, n_x, n_y, max_size = 0;
+ int *clust_index, *clust_size, max_clust_size, max_clust_ind, nav, nhisto;
+ t_rgb rlo = { 1.0, 1.0, 1.0 };
+
+ clear_trxframe(&fr, TRUE);
+ sprintf(timebuf, "Time (%s)", output_env_get_time_unit(oenv));
+ tf = output_env_get_time_factor(oenv);
+ fp = xvgropen(ncl, "Number of clusters", timebuf, "N", oenv);
+ gp = xvgropen(acl, "Average cluster size", timebuf, "#molecules", oenv);
+ hp = xvgropen(mcl, "Max cluster size", timebuf, "#molecules", oenv);
+ tp = xvgropen(tempf, "Temperature of largest cluster", timebuf, "T (K)",
+ oenv);
+
+ if (!read_first_frame(oenv, &status, trx, &fr, TRX_NEED_X | TRX_READ_V))
+ {
+ gmx_file(trx);
+ }
+
+ natoms = fr.natoms;
+ x = fr.x;
+
+ if (tpr)
+ {
+ snew(mtop, 1);
+ read_tpxheader(tpr, &tpxh, TRUE, &version, &generation);
+ if (tpxh.natoms != natoms)
+ {
+ gmx_fatal(FARGS, "tpr (%d atoms) and xtc (%d atoms) do not match!",
+ tpxh.natoms, natoms);
+ }
+ ePBC = read_tpx(tpr, NULL, NULL, &natoms, NULL, NULL, NULL, mtop);
+ }
+ if (ndf <= -1)
+ {
+ tfac = 1;
}
- /* Analyse velocities, if present */
- if (fr.bV) {
- if (!tpr) {
- if (bTPRwarn) {
- printf("You need a [TT].tpr[tt] file to analyse temperatures\n");
- bTPRwarn = FALSE;
- }
- }
- else {
- v = fr.v;
- /* Loop over clusters and for each cluster compute 1/2 m v^2 */
- if (max_clust_ind >= 0) {
- ekin = 0;
- for(i=0; (i<nindex); i++)
- if (clust_index[i] == max_clust_ind) {
- ai = index[i];
- gmx_mtop_atomnr_to_atom(alook,ai,&atom);
- ekin += 0.5*atom->m*iprod(v[ai],v[ai]);
- }
- temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
- fprintf(tp,"%10.3f %10.3f\n",fr.time,temp);
- }
- }
+ else
+ {
+ tfac = ndf/(3.0*natoms);
}
- nframe++;
- } while (read_next_frame(oenv,status,&fr));
- close_trx(status);
- ffclose(fp);
- ffclose(gp);
- ffclose(hp);
- ffclose(tp);
-
- gmx_mtop_atomlookup_destroy(alook);
-
- if (max_clust_ind >= 0) {
- fp = ffopen(mcn,"w");
- fprintf(fp,"[ max_clust ]\n");
- for(i=0; (i<nindex); i++)
- if (clust_index[i] == max_clust_ind) {
- if (bMol) {
- for(j=mols->index[i]; (j<mols->index[i+1]); j++)
- fprintf(fp,"%d\n",j+1);
- }
- else {
- fprintf(fp,"%d\n",index[i]+1);
- }
- }
+
+ if (bMol)
+ {
+ if (ndx)
+ {
+ printf("Using molecules rather than atoms. Not reading index file %s\n",
+ ndx);
+ }
+ mols = &(mtop->mols);
+
+ /* Make dummy index */
+ nindex = mols->nr;
+ snew(index, nindex);
+ for (i = 0; (i < nindex); i++)
+ {
+ index[i] = i;
+ }
+ gname = strdup("mols");
+ }
+ else
+ {
+ rd_index(ndx, 1, &nindex, &index, &gname);
+ }
+
+ alook = gmx_mtop_atomlookup_init(mtop);
+
+ snew(clust_index, nindex);
+ snew(clust_size, nindex);
+ cut2 = cut*cut;
+ nframe = 0;
+ n_x = 0;
+ snew(t_y, nindex);
+ for (i = 0; (i < nindex); i++)
+ {
+ t_y[i] = i+1;
+ }
+ max_clust_size = 1;
+ max_clust_ind = -1;
+ do
+ {
+ if ((nskip == 0) || ((nskip > 0) && ((nframe % nskip) == 0)))
+ {
+ if (bPBC)
+ {
+ set_pbc(&pbc, ePBC, fr.box);
+ }
+ max_clust_size = 1;
+ max_clust_ind = -1;
+
+ /* Put all atoms/molecules in their own cluster, with size 1 */
+ for (i = 0; (i < nindex); i++)
+ {
+ /* Cluster index is indexed with atom index number */
+ clust_index[i] = i;
+ /* Cluster size is indexed with cluster number */
+ clust_size[i] = 1;
+ }
+
+ /* Loop over atoms */
+ for (i = 0; (i < nindex); i++)
+ {
+ ai = index[i];
+ ci = clust_index[i];
+
+ /* Loop over atoms (only half a matrix) */
+ for (j = i+1; (j < nindex); j++)
+ {
+ cj = clust_index[j];
+
+ /* If they are not in the same cluster already */
+ if (ci != cj)
+ {
+ aj = index[j];
+
+ /* Compute distance */
+ if (bMol)
+ {
+ bSame = FALSE;
+ for (ii = mols->index[ai]; !bSame && (ii < mols->index[ai+1]); ii++)
+ {
+ for (jj = mols->index[aj]; !bSame && (jj < mols->index[aj+1]); jj++)
+ {
+ if (bPBC)
+ {
+ pbc_dx(&pbc, x[ii], x[jj], dx);
+ }
+ else
+ {
+ rvec_sub(x[ii], x[jj], dx);
+ }
+ dx2 = iprod(dx, dx);
+ bSame = (dx2 < cut2);
+ }
+ }
+ }
+ else
+ {
+ if (bPBC)
+ {
+ pbc_dx(&pbc, x[ai], x[aj], dx);
+ }
+ else
+ {
+ rvec_sub(x[ai], x[aj], dx);
+ }
+ dx2 = iprod(dx, dx);
+ bSame = (dx2 < cut2);
+ }
+ /* If distance less than cut-off */
+ if (bSame)
+ {
+ /* Merge clusters: check for all atoms whether they are in
+ * cluster cj and if so, put them in ci
+ */
+ for (k = 0; (k < nindex); k++)
+ {
+ if (clust_index[k] == cj)
+ {
+ if (clust_size[cj] <= 0)
+ {
+ gmx_fatal(FARGS, "negative cluster size %d for element %d",
+ clust_size[cj], cj);
+ }
+ clust_size[cj]--;
+ clust_index[k] = ci;
+ clust_size[ci]++;
+ }
+ }
+ }
+ }
+ }
+ }
+ n_x++;
+ srenew(t_x, n_x);
+ t_x[n_x-1] = fr.time*tf;
+ srenew(cs_dist, n_x);
+ snew(cs_dist[n_x-1], nindex);
+ nclust = 0;
+ cav = 0;
+ nav = 0;
+ for (i = 0; (i < nindex); i++)
+ {
+ ci = clust_size[i];
+ if (ci > max_clust_size)
+ {
+ max_clust_size = ci;
+ max_clust_ind = i;
+ }
+ if (ci > 0)
+ {
+ nclust++;
+ cs_dist[n_x-1][ci-1] += 1.0;
+ max_size = max(max_size, ci);
+ if (ci > 1)
+ {
+ cav += ci;
+ nav++;
+ }
+ }
+ }
+ fprintf(fp, "%14.6e %10d\n", fr.time, nclust);
+ if (nav > 0)
+ {
+ fprintf(gp, "%14.6e %10.3f\n", fr.time, cav/nav);
+ }
+ fprintf(hp, "%14.6e %10d\n", fr.time, max_clust_size);
+ }
+ /* Analyse velocities, if present */
+ if (fr.bV)
+ {
+ if (!tpr)
+ {
+ if (bTPRwarn)
+ {
+ printf("You need a [TT].tpr[tt] file to analyse temperatures\n");
+ bTPRwarn = FALSE;
+ }
+ }
+ else
+ {
+ v = fr.v;
+ /* Loop over clusters and for each cluster compute 1/2 m v^2 */
+ if (max_clust_ind >= 0)
+ {
+ ekin = 0;
+ for (i = 0; (i < nindex); i++)
+ {
+ if (clust_index[i] == max_clust_ind)
+ {
+ ai = index[i];
+ gmx_mtop_atomnr_to_atom(alook, ai, &atom);
+ ekin += 0.5*atom->m*iprod(v[ai], v[ai]);
+ }
+ }
+ temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
+ fprintf(tp, "%10.3f %10.3f\n", fr.time, temp);
+ }
+ }
+ }
+ nframe++;
+ }
+ while (read_next_frame(oenv, status, &fr));
+ close_trx(status);
ffclose(fp);
- }
-
- /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
- fp = xvgropen(histo,"Cluster size distribution","Cluster size","()",oenv);
- nhisto = 0;
- fprintf(fp,"%5d %8.3f\n",0,0.0);
- for(j=0; (j<max_size); j++) {
- real nelem = 0;
- for(i=0; (i<n_x); i++)
- nelem += cs_dist[i][j];
- fprintf(fp,"%5d %8.3f\n",j+1,nelem/n_x);
- nhisto += (int)((j+1)*nelem/n_x);
- }
- fprintf(fp,"%5d %8.3f\n",j+1,0.0);
- ffclose(fp);
-
- fprintf(stderr,"Total number of atoms in clusters = %d\n",nhisto);
-
- /* Look for the smallest entry that is not zero
- * This will make that zero is white, and not zero is coloured.
- */
- cmid = 100.0;
- cmax = 0.0;
- for(i=0; (i<n_x); i++)
- for(j=0; (j<max_size); j++) {
- if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
- cmid = cs_dist[i][j];
- cmax = max(cs_dist[i][j],cmax);
+ ffclose(gp);
+ ffclose(hp);
+ ffclose(tp);
+
+ gmx_mtop_atomlookup_destroy(alook);
+
+ if (max_clust_ind >= 0)
+ {
+ fp = ffopen(mcn, "w");
+ fprintf(fp, "[ max_clust ]\n");
+ for (i = 0; (i < nindex); i++)
+ {
+ if (clust_index[i] == max_clust_ind)
+ {
+ if (bMol)
+ {
+ for (j = mols->index[i]; (j < mols->index[i+1]); j++)
+ {
+ fprintf(fp, "%d\n", j+1);
+ }
+ }
+ else
+ {
+ fprintf(fp, "%d\n", index[i]+1);
+ }
+ }
+ }
+ ffclose(fp);
}
- fprintf(stderr,"cmid: %g, cmax: %g, max_size: %d\n",cmid,cmax,max_size);
- cmid = 1;
- fp = ffopen(xpm,"w");
- write_xpm3(fp,0,"Cluster size distribution","# clusters",timebuf,"Size",
- n_x,max_size,t_x,t_y,cs_dist,0,cmid,cmax,
- rlo,rmid,rhi,&nlevels);
- ffclose(fp);
- cmid = 100.0;
- cmax = 0.0;
- for(i=0; (i<n_x); i++)
- for(j=0; (j<max_size); j++) {
- cs_dist[i][j] *= (j+1);
- if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
- cmid = cs_dist[i][j];
- cmax = max(cs_dist[i][j],cmax);
+
+ /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
+ fp = xvgropen(histo, "Cluster size distribution", "Cluster size", "()", oenv);
+ nhisto = 0;
+ fprintf(fp, "%5d %8.3f\n", 0, 0.0);
+ for (j = 0; (j < max_size); j++)
+ {
+ real nelem = 0;
+ for (i = 0; (i < n_x); i++)
+ {
+ nelem += cs_dist[i][j];
+ }
+ fprintf(fp, "%5d %8.3f\n", j+1, nelem/n_x);
+ nhisto += (int)((j+1)*nelem/n_x);
+ }
+ fprintf(fp, "%5d %8.3f\n", j+1, 0.0);
+ ffclose(fp);
+
+ fprintf(stderr, "Total number of atoms in clusters = %d\n", nhisto);
+
+ /* Look for the smallest entry that is not zero
+ * This will make that zero is white, and not zero is coloured.
+ */
+ cmid = 100.0;
+ cmax = 0.0;
+ for (i = 0; (i < n_x); i++)
+ {
+ for (j = 0; (j < max_size); j++)
+ {
+ if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
+ {
+ cmid = cs_dist[i][j];
+ }
+ cmax = max(cs_dist[i][j], cmax);
+ }
}
- fprintf(stderr,"cmid: %g, cmax: %g, max_size: %d\n",cmid,cmax,max_size);
- fp = ffopen(xpmw,"w");
- write_xpm3(fp,0,"Weighted cluster size distribution","Fraction",timebuf,
- "Size", n_x,max_size,t_x,t_y,cs_dist,0,cmid,cmax,
- rlo,rmid,rhi,&nlevels);
- ffclose(fp);
-
- sfree(clust_index);
- sfree(clust_size);
- sfree(index);
+ fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
+ cmid = 1;
+ fp = ffopen(xpm, "w");
+ write_xpm3(fp, 0, "Cluster size distribution", "# clusters", timebuf, "Size",
+ n_x, max_size, t_x, t_y, cs_dist, 0, cmid, cmax,
+ rlo, rmid, rhi, &nlevels);
+ ffclose(fp);
+ cmid = 100.0;
+ cmax = 0.0;
+ for (i = 0; (i < n_x); i++)
+ {
+ for (j = 0; (j < max_size); j++)
+ {
+ cs_dist[i][j] *= (j+1);
+ if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
+ {
+ cmid = cs_dist[i][j];
+ }
+ cmax = max(cs_dist[i][j], cmax);
+ }
+ }
+ fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
+ fp = ffopen(xpmw, "w");
+ write_xpm3(fp, 0, "Weighted cluster size distribution", "Fraction", timebuf,
+ "Size", n_x, max_size, t_x, t_y, cs_dist, 0, cmid, cmax,
+ rlo, rmid, rhi, &nlevels);
+ ffclose(fp);
+
+ sfree(clust_index);
+ sfree(clust_size);
+ sfree(index);
}
-int gmx_clustsize(int argc,char *argv[])
+int gmx_clustsize(int argc, char *argv[])
{
- const char *desc[] = {
- "This program computes the size distributions of molecular/atomic clusters in",
- "the gas phase. The output is given in the form of an [TT].xpm[tt] file.",
- "The total number of clusters is written to an [TT].xvg[tt] file.[PAR]",
- "When the [TT]-mol[tt] option is given clusters will be made out of",
- "molecules rather than atoms, which allows clustering of large molecules.",
- "In this case an index file would still contain atom numbers",
- "or your calculation will die with a SEGV.[PAR]",
- "When velocities are present in your trajectory, the temperature of",
- "the largest cluster will be printed in a separate [TT].xvg[tt] file assuming",
- "that the particles are free to move. If you are using constraints,",
- "please correct the temperature. For instance water simulated with SHAKE",
- "or SETTLE will yield a temperature that is 1.5 times too low. You can",
- "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal",
- "of center of mass motion into account.[PAR]",
- "The [TT]-mc[tt] option will produce an index file containing the",
- "atom numbers of the largest cluster."
- };
-
- static real cutoff = 0.35;
- static int nskip = 0;
- static int nlevels = 20;
- static int ndf = -1;
- static gmx_bool bMol = FALSE;
- static gmx_bool bPBC = TRUE;
- static rvec rlo = { 1.0, 1.0, 0.0 };
- static rvec rhi = { 0.0, 0.0, 1.0 };
-
- output_env_t oenv;
-
- t_pargs pa[] = {
- { "-cut", FALSE, etREAL, {&cutoff},
- "Largest distance (nm) to be considered in a cluster" },
- { "-mol", FALSE, etBOOL, {&bMol},
- "Cluster molecules rather than atoms (needs [TT].tpr[tt] file)" },
- { "-pbc", FALSE, etBOOL, {&bPBC},
- "Use periodic boundary conditions" },
- { "-nskip", FALSE, etINT, {&nskip},
- "Number of frames to skip between writing" },
- { "-nlevels", FALSE, etINT, {&nlevels},
- "Number of levels of grey in [TT].xpm[tt] output" },
- { "-ndf", FALSE, etINT, {&ndf},
- "Number of degrees of freedom of the entire system for temperature calculation. If not set, the number of atoms times three is used." },
- { "-rgblo", FALSE, etRVEC, {rlo},
- "RGB values for the color of the lowest occupied cluster size" },
- { "-rgbhi", FALSE, etRVEC, {rhi},
- "RGB values for the color of the highest occupied cluster size" }
- };
+ const char *desc[] = {
+ "This program computes the size distributions of molecular/atomic clusters in",
+ "the gas phase. The output is given in the form of an [TT].xpm[tt] file.",
+ "The total number of clusters is written to an [TT].xvg[tt] file.[PAR]",
+ "When the [TT]-mol[tt] option is given clusters will be made out of",
+ "molecules rather than atoms, which allows clustering of large molecules.",
+ "In this case an index file would still contain atom numbers",
+ "or your calculation will die with a SEGV.[PAR]",
+ "When velocities are present in your trajectory, the temperature of",
+ "the largest cluster will be printed in a separate [TT].xvg[tt] file assuming",
+ "that the particles are free to move. If you are using constraints,",
+ "please correct the temperature. For instance water simulated with SHAKE",
+ "or SETTLE will yield a temperature that is 1.5 times too low. You can",
+ "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal",
+ "of center of mass motion into account.[PAR]",
+ "The [TT]-mc[tt] option will produce an index file containing the",
+ "atom numbers of the largest cluster."
+ };
+
+ static real cutoff = 0.35;
+ static int nskip = 0;
+ static int nlevels = 20;
+ static int ndf = -1;
+ static gmx_bool bMol = FALSE;
+ static gmx_bool bPBC = TRUE;
+ static rvec rlo = { 1.0, 1.0, 0.0 };
+ static rvec rhi = { 0.0, 0.0, 1.0 };
+
+ output_env_t oenv;
+
+ t_pargs pa[] = {
+ { "-cut", FALSE, etREAL, {&cutoff},
+ "Largest distance (nm) to be considered in a cluster" },
+ { "-mol", FALSE, etBOOL, {&bMol},
+ "Cluster molecules rather than atoms (needs [TT].tpr[tt] file)" },
+ { "-pbc", FALSE, etBOOL, {&bPBC},
+ "Use periodic boundary conditions" },
+ { "-nskip", FALSE, etINT, {&nskip},
+ "Number of frames to skip between writing" },
+ { "-nlevels", FALSE, etINT, {&nlevels},
+ "Number of levels of grey in [TT].xpm[tt] output" },
+ { "-ndf", FALSE, etINT, {&ndf},
+ "Number of degrees of freedom of the entire system for temperature calculation. If not set, the number of atoms times three is used." },
+ { "-rgblo", FALSE, etRVEC, {rlo},
+ "RGB values for the color of the lowest occupied cluster size" },
+ { "-rgbhi", FALSE, etRVEC, {rhi},
+ "RGB values for the color of the highest occupied cluster size" }
+ };
#define NPA asize(pa)
- const char *fnNDX,*fnTPR;
- gmx_bool bSQ,bRDF;
- t_rgb rgblo,rgbhi;
-
- t_filenm fnm[] = {
- { efTRX, "-f", NULL, ffREAD },
- { efTPR, NULL, NULL, ffOPTRD },
- { efNDX, NULL, NULL, ffOPTRD },
- { efXPM, "-o", "csize", ffWRITE },
- { efXPM, "-ow","csizew", ffWRITE },
- { efXVG, "-nc","nclust", ffWRITE },
- { efXVG, "-mc","maxclust", ffWRITE },
- { efXVG, "-ac","avclust", ffWRITE },
- { efXVG, "-hc","histo-clust", ffWRITE },
- { efXVG, "-temp","temp", ffOPTWR },
- { efNDX, "-mcn", "maxclust", ffOPTWR }
- };
+ const char *fnNDX, *fnTPR;
+ gmx_bool bSQ, bRDF;
+ t_rgb rgblo, rgbhi;
+
+ t_filenm fnm[] = {
+ { efTRX, "-f", NULL, ffREAD },
+ { efTPR, NULL, NULL, ffOPTRD },
+ { efNDX, NULL, NULL, ffOPTRD },
+ { efXPM, "-o", "csize", ffWRITE },
+ { efXPM, "-ow", "csizew", ffWRITE },
+ { efXVG, "-nc", "nclust", ffWRITE },
+ { efXVG, "-mc", "maxclust", ffWRITE },
+ { efXVG, "-ac", "avclust", ffWRITE },
+ { efXVG, "-hc", "histo-clust", ffWRITE },
+ { efXVG, "-temp", "temp", ffOPTWR },
+ { efNDX, "-mcn", "maxclust", ffOPTWR }
+ };
#define NFILE asize(fnm)
-
- parse_common_args(&argc,argv,
- PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
- NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
-
- fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
- rgblo.r = rlo[XX],rgblo.g = rlo[YY],rgblo.b = rlo[ZZ];
- rgbhi.r = rhi[XX],rgbhi.g = rhi[YY],rgbhi.b = rhi[ZZ];
-
- fnTPR = ftp2fn_null(efTPR,NFILE,fnm);
- if (bMol && !fnTPR)
- gmx_fatal(FARGS,"You need a tpr file for the -mol option");
-
- clust_size(fnNDX,ftp2fn(efTRX,NFILE,fnm),opt2fn("-o",NFILE,fnm),
- opt2fn("-ow",NFILE,fnm),
- opt2fn("-nc",NFILE,fnm),opt2fn("-ac",NFILE,fnm),
- opt2fn("-mc",NFILE,fnm),opt2fn("-hc",NFILE,fnm),
- opt2fn("-temp",NFILE,fnm),opt2fn("-mcn",NFILE,fnm),
- bMol,bPBC,fnTPR,
- cutoff,nskip,nlevels,rgblo,rgbhi,ndf,oenv);
-
- thanx(stderr);
-
- return 0;
+
+ parse_common_args(&argc, argv,
+ PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
+ NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
+
+ fnNDX = ftp2fn_null(efNDX, NFILE, fnm);
+ rgblo.r = rlo[XX], rgblo.g = rlo[YY], rgblo.b = rlo[ZZ];
+ rgbhi.r = rhi[XX], rgbhi.g = rhi[YY], rgbhi.b = rhi[ZZ];
+
+ fnTPR = ftp2fn_null(efTPR, NFILE, fnm);
+ if (bMol && !fnTPR)
+ {
+ gmx_fatal(FARGS, "You need a tpr file for the -mol option");
+ }
+
+ clust_size(fnNDX, ftp2fn(efTRX, NFILE, fnm), opt2fn("-o", NFILE, fnm),
+ opt2fn("-ow", NFILE, fnm),
+ opt2fn("-nc", NFILE, fnm), opt2fn("-ac", NFILE, fnm),
+ opt2fn("-mc", NFILE, fnm), opt2fn("-hc", NFILE, fnm),
+ opt2fn("-temp", NFILE, fnm), opt2fn("-mcn", NFILE, fnm),
+ bMol, bPBC, fnTPR,
+ cutoff, nskip, nlevels, rgblo, rgbhi, ndf, oenv);
+
+ thanx(stderr);
+
+ return 0;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff