-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*- */
/*
+ * This file is part of the GROMACS molecular simulation package.
*
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
* 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.
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
*
- * 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.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
- * For more info, check our website at http://www.gromacs.org
+ * If you want to redistribute modifications to GROMACS, 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 http://www.gromacs.org.
*
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
/*! \internal \file
* \author Jochen Hub <jhub@gwdg.de>
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "config.h"
+
+#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
+#include <string.h>
+
#include <sstream>
-#include "statutil.h"
-#include "typedefs.h"
-#include "smalloc.h"
-#include "vec.h"
-#include "copyrite.h"
-#include "statutil.h"
-#include "tpxio.h"
-#include "names.h"
-#include "gmx_random.h"
-#include "gmx_ana.h"
-#include "macros.h"
-
-#include "string2.h"
-#include "xvgr.h"
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/random/random.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
//! longest file names allowed in input files
#define WHAM_MAXFILELEN 2048
+/*! \brief
+ * x-axis legend for output files
+ */
+static const char *xlabel = "\\xx\\f{} (nm)";
+
/*! \brief
* enum for energy units
*/
* \name Using umbrella pull code since gromacs 4.x
*/
/*!\{*/
- int npullgrps; //!< nr of pull groups in tpr file
- int pull_geometry; //!< such as distance, position
- ivec pull_dim; //!< pull dimension with geometry distance
- int pull_ndim; //!< nr of pull_dim != 0
- real *k; //!< force constants in tpr file
- rvec *init_dist; //!< reference displacements
- real *umbInitDist; //!< reference displacement in umbrella direction
+ int npullcrds; //!< nr of pull coordinates in tpr file
+ int pull_geometry; //!< such as distance, direction
+ ivec pull_dim; //!< pull dimension with geometry distance
+ int pull_ndim; //!< nr of pull_dim != 0
+ gmx_bool bPrintRef; //!< Coordinates of reference groups written to pullx.xvg ?
+ real *k; //!< force constants in tpr file
+ real *init_dist; //!< reference displacements
+ real *umbInitDist; //!< reference displacement in umbrella direction
/*!\}*/
/*!
* \name Using PDO files common until gromacs 3.x
snew(fn, strlen(fnhist)+10);
snew(buf, strlen(fnhist)+10);
sprintf(fn, "%s_cumul.xvg", strncpy(buf, fnhist, strlen(fnhist)-4));
- fp = xvgropen(fn, "CDFs of umbrella windows", "z", "CDF", opt->oenv);
+ fp = xvgropen(fn, "CDFs of umbrella windows", xlabel, "CDF", opt->oenv);
}
nbin = opt->bins;
fprintf(fp, "\n");
}
printf("Wrote cumulative distribution functions to %s\n", fn);
- ffclose(fp);
+ gmx_ffclose(fp);
sfree(fn);
sfree(buf);
}
}
else
{
- fn = strdup(fnhist);
+ fn = gmx_strdup(fnhist);
strcpy(title, "Umbrella histograms");
}
- fp = xvgropen(fn, title, "z", "count", opt->oenv);
+ fp = xvgropen(fn, title, xlabel, "count", opt->oenv);
bins = opt->bins;
/* Write histograms */
fprintf(fp, "\n");
}
- ffclose(fp);
+ gmx_ffclose(fp);
printf("Wrote %s\n", fn);
if (bs_index >= 0)
{
}
/* do bootstrapping */
- fp = xvgropen(fnprof, "Boot strap profiles", "z", ylabel, opt->oenv);
+ fp = xvgropen(fnprof, "Boot strap profiles", xlabel, ylabel, opt->oenv);
for (ib = 0; ib < opt->nBootStrap; ib++)
{
printf(" *******************************************\n"
bsProfiles_av2[i] += tmp*tmp;
fprintf(fp, "%e\t%e\n", (i+0.5)*opt->dz+opt->min, tmp);
}
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
}
- ffclose(fp);
+ gmx_ffclose(fp);
/* write average and stddev */
- fp = xvgropen(fnres, "Average and stddev from bootstrapping", "z", ylabel, opt->oenv);
- fprintf(fp, "@TYPE xydy\n");
+ fp = xvgropen(fnres, "Average and stddev from bootstrapping", xlabel, ylabel, opt->oenv);
+ if (output_env_get_print_xvgr_codes(opt->oenv))
+ {
+ fprintf(fp, "@TYPE xydy\n");
+ }
for (i = 0; i < opt->bins; i++)
{
bsProfiles_av [i] /= opt->nBootStrap;
stddev = (tmp >= 0.) ? sqrt(tmp) : 0.; /* Catch rouding errors */
fprintf(fp, "%e\t%e\t%e\n", (i+0.5)*opt->dz+opt->min, bsProfiles_av [i], stddev);
}
- ffclose(fp);
+ gmx_ffclose(fp);
printf("Wrote boot strap result to %s\n", fnres);
}
int nread, sizenow, i, block = 1;
FILE *fp;
- fp = ffopen(fn, "r");
+ fp = gmx_ffopen(fn, "r");
nread = 0;
sizenow = 0;
while (fgets(tmp, sizeof(tmp), fp) != NULL)
}
else
{
- pipe = ffopen(fn, "r");
+ pipe = gmx_ffopen(fn, "r");
*bPipeOpen = FALSE;
}
#ifdef HAVE_PIPES
pclose(fp);
#else
- ffclose(fp);
+ gmx_ffclose(fp);
#endif
}
}
else
{
- ffclose(file);
+ gmx_ffclose(file);
}
}
printf("\n");
}
else
{
- ffclose(file);
+ gmx_ffclose(file);
}
}
printf("\n");
void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions *opt)
{
t_inputrec ir;
- int i, ngrp, d;
+ int i, ncrd;
t_state state;
static int first = 1;
}
/* nr of pull groups */
- ngrp = ir.pull->ngrp;
- if (ngrp < 1)
+ ncrd = ir.pull->ncoord;
+ if (ncrd < 1)
{
- gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull groups\n", ngrp);
+ gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull coordinates\n", ncrd);
}
- header->npullgrps = ir.pull->ngrp;
+ header->npullcrds = ir.pull->ncoord;
header->pull_geometry = ir.pull->eGeom;
+ header->bPrintRef = ir.pull->bPrintRef;
copy_ivec(ir.pull->dim, header->pull_dim);
header->pull_ndim = header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2];
- if (header->pull_geometry == epullgPOS && header->pull_ndim > 1)
- {
- gmx_fatal(FARGS, "Found pull geometry 'position' and more than 1 pull dimension (%d).\n"
- "Hence, the pull potential does not correspond to a one-dimensional umbrella potential.\n"
- "If you have some special umbrella setup you may want to write your own pdo files\n"
- "and feed them into g_wham. Check g_wham -h !\n", header->pull_ndim);
- }
- snew(header->k, ngrp);
- snew(header->init_dist, ngrp);
- snew(header->umbInitDist, ngrp);
+ snew(header->k, ncrd);
+ snew(header->init_dist, ncrd);
+ snew(header->umbInitDist, ncrd);
/* only z-direction with epullgCYL? */
if (header->pull_geometry == epullgCYL)
}
}
- for (i = 0; i < ngrp; i++)
+ for (i = 0; i < ncrd; i++)
{
- header->k[i] = ir.pull->grp[i+1].k;
+ header->k[i] = ir.pull->coord[i].k;
if (header->k[i] == 0.0)
{
- gmx_fatal(FARGS, "Pull group %d has force constant of of 0.0 in %s.\n"
+ gmx_fatal(FARGS, "Pull coordinate %d has force constant of of 0.0 in %s.\n"
"That doesn't seem to be an Umbrella tpr.\n",
i, fn);
}
- copy_rvec(ir.pull->grp[i+1].init, header->init_dist[i]);
+ header->init_dist[i] = ir.pull->coord[i].init;
/* initial distance to reference */
switch (header->pull_geometry)
{
- case epullgPOS:
- for (d = 0; d < DIM; d++)
- {
- if (header->pull_dim[d])
- {
- header->umbInitDist[i] = header->init_dist[i][d];
- }
- }
- break;
case epullgCYL:
- /* umbrella distance stored in init_dist[i][0] for geometry cylinder (not in ...[i][ZZ]) */
+ /* umbrella distance stored in init_dist[i] for geometry cylinder (not in ...[i][ZZ]) */
case epullgDIST:
case epullgDIR:
case epullgDIRPBC:
- header->umbInitDist[i] = header->init_dist[i][0];
+ header->umbInitDist[i] = header->init_dist[i];
break;
default:
gmx_fatal(FARGS, "Pull geometry %s not supported\n", epullg_names[header->pull_geometry]);
if (opt->verbose || first)
{
- printf("File %s, %d groups, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n",
- fn, header->npullgrps, epullg_names[header->pull_geometry],
+ printf("File %s, %d coordinates, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n"
+ "\tPull group coordinates%s expected in pullx files.\n",
+ fn, header->npullcrds, epullg_names[header->pull_geometry],
int2YN(header->pull_dim[0]), int2YN(header->pull_dim[1]), int2YN(header->pull_dim[2]),
- header->pull_ndim);
- for (i = 0; i < ngrp; i++)
+ header->pull_ndim, (header->bPrintRef ? "" : " not"));
+ for (i = 0; i < ncrd; i++)
{
- printf("\tgrp %d) k = %-5g position = %g\n", i, header->k[i], header->umbInitDist[i]);
+ printf("\tcrd %d) k = %-5g position = %g\n", i, header->k[i], header->umbInitDist[i]);
}
}
if (!opt->verbose && first)
{
- printf("\tUse option -v to see this output for all input tpr files\n");
+ printf("\tUse option -v to see this output for all input tpr files\n\n");
}
first = 0;
t_groupselection *groupsel)
{
double **y = 0, pos = 0., t, force, time0 = 0., dt;
- int ny, nt, bins, ibin, i, g, gUsed, dstep = 1, nColPerGrp, nColRefOnce, nColRefEachGrp, nColExpect, ntot;
+ int ny, nt, bins, ibin, i, g, gUsed, dstep = 1, nColPerCrd, nColRefEachCrd, nColExpect, ntot, column;
real min, max, minfound = 1e20, maxfound = -1e20;
gmx_bool dt_ok, timeok, bHaveForce;
const char *quantity;
static gmx_bool bFirst = TRUE;
/*
- in force output pullf.xvg:
- No reference, one column per pull group
- in position output pullx.xvg (not cylinder)
- ndim reference, ndim columns per pull group
- in position output pullx.xvg (in geometry cylinder):
- ndim*2 columns per pull group (ndim for ref, ndim for group)
+ * Data columns in pull output:
+ * - in force output pullf.xvg:
+ * No reference columns, one column per pull coordinate
+ *
+ * - in position output pullx.xvg
+ * bPrintRef == TRUE: for each pull coordinate: ndim reference columns, and ndim dx columns
+ * bPrintRef == FALSE: for each pull coordinate: no reference columns, but ndim dx columns
*/
- nColPerGrp = opt->bPullx ? header->pull_ndim : 1;
+ nColPerCrd = opt->bPullx ? header->pull_ndim : 1;
quantity = opt->bPullx ? "position" : "force";
- if (opt->bPullx)
+ if (opt->bPullx && header->bPrintRef)
{
- if (header->pull_geometry == epullgCYL)
- {
- /* Geometry cylinder -> reference group before each pull group */
- nColRefEachGrp = header->pull_ndim;
- nColRefOnce = 0;
- }
- else
- {
- /* Geometry NOT cylinder -> reference group only once after time column */
- nColRefEachGrp = 0;
- nColRefOnce = header->pull_ndim;
- }
+ nColRefEachCrd = header->pull_ndim;
}
- else /* read forces, no reference groups */
+ else
{
- nColRefEachGrp = 0;
- nColRefOnce = 0;
+ nColRefEachCrd = 0;
}
- nColExpect = 1 + nColRefOnce + header->npullgrps*(nColRefEachGrp+nColPerGrp);
+ nColExpect = 1 + header->npullcrds*(nColRefEachCrd+nColPerCrd);
bHaveForce = opt->bPullf;
/* With geometry "distance" or "distance_periodic", only force reading is supported so far.
- That avoids the somewhat tedious extraction of the right columns from the pullx files
- to compute the distances projection on the vector. Sorry for the laziness. */
+ That avoids the somewhat tedious extraction of the reaction coordinate from the pullx files.
+ Sorry for the laziness, this is a To-do. */
if ( (header->pull_geometry == epullgDIR || header->pull_geometry == epullgDIRPBC)
&& opt->bPullx)
{
bHaveForce ? "force" : "position", epullg_names[header->pull_geometry],
header->pull_ndim);
printf("Expecting these columns in pull file:\n"
- "\t%d reference columns for all pull groups together\n"
- "\t%d reference columns for each individual pull group\n"
- "\t%d data columns for each pull group\n", nColRefOnce, nColRefEachGrp, nColPerGrp);
- printf("With %d pull groups, expect %d columns (including the time column)\n", header->npullgrps, nColExpect);
+ "\t%d reference columns for each individual pull coordinate\n"
+ "\t%d data columns for each pull coordinate\n", nColRefEachCrd, nColPerCrd);
+ printf("With %d pull groups, expect %d columns (including the time column)\n", header->npullcrds, nColExpect);
bFirst = FALSE;
}
if (ny != nColExpect)
{
- gmx_fatal(FARGS, "Found %d pull groups in %s,\n but %d data columns in %s (expected %d)\n"
+ gmx_fatal(FARGS, "Found %d pull coodinates in %s,\n but %d data columns in %s (expected %d)\n"
"\nMaybe you confused options -ix and -if ?\n",
- header->npullgrps, fntpr, ny-1, fn, nColExpect-1);
+ header->npullcrds, fntpr, ny-1, fn, nColExpect-1);
}
if (opt->verbose)
{
- printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerGrp, quantity, fn);
+ printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerCrd, quantity, fn);
}
if (!bGetMinMax)
if (groupsel)
{
/* Use only groups selected with option -is file */
- if (header->npullgrps != groupsel->n)
+ if (header->npullcrds != groupsel->n)
{
gmx_fatal(FARGS, "tpr file contains %d pull groups, but expected %d from group selection file\n",
- header->npullgrps, groupsel->n);
+ header->npullcrds, groupsel->n);
}
window->nPull = groupsel->nUse;
}
else
{
- window->nPull = header->npullgrps;
+ window->nPull = header->npullcrds;
}
window->nBin = bins;
/* Copying umbrella center and force const is more involved since not
all pull groups from header (tpr file) may be used in window variable */
- for (g = 0, gUsed = 0; g < header->npullgrps; ++g)
+ for (g = 0, gUsed = 0; g < header->npullcrds; ++g)
{
if (groupsel && (groupsel->bUse[g] == FALSE))
{
/*if (opt->verbose)
printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
t,opt->tmin, opt->tmax, dt_ok,timeok); */
-
if (timeok)
{
/* Note: if groupsel == NULL:
* only groups with groupsel.bUse[g]==TRUE are stored. gUsed is not always equal g
*/
gUsed = -1;
- for (g = 0; g < header->npullgrps; ++g)
+ for (g = 0; g < header->npullcrds; ++g)
{
/* was this group selected for application in WHAM? */
if (groupsel && (groupsel->bUse[g] == FALSE))
if (bHaveForce)
{
- /* y has 1 time column y[0] and one column per force y[1],...,y[nGrps] */
+ /* y has 1 time column y[0] and one column per force y[1],...,y[nCrds] */
force = y[g+1][i];
pos = -force/header->k[g] + header->umbInitDist[g];
}
else
{
+ /* pick the right column from:
+ * time ref1[ndim] coord1[ndim] ref2[ndim] coord2[ndim] ...
+ */
+ column = 1 + nColRefEachCrd + g*(nColRefEachCrd+nColPerCrd);
switch (header->pull_geometry)
{
case epullgDIST:
- /* y has 1 time column y[0] and nColPerGrps columns per pull group;
- Distance to reference: */
- /* pos=dist_ndim(y+1+nColRef+g*nColPerGrp,header->pull_ndim,i); gmx 4.0 */
- pos = dist_ndim(y + 1 + nColRefOnce + g*nColPerGrp, header->pull_ndim, i);
+ pos = dist_ndim(y + column, header->pull_ndim, i);
break;
- case epullgPOS:
- /* Columns
- Time ref[ndim] group1[ndim] group2[ndim] ... */
case epullgCYL:
- /* Columns
- Time ref1[ndim] group1[ndim] ref2[ndim] group2[ndim] ... */
-
- /* * with geometry==position, we have the reference once (nColRefOnce==ndim), but
- no extra reference group columns before each group (nColRefEachGrp==0)
-
- * with geometry==cylinder, we have no initial ref group column (nColRefOnce==0),
- but ndim ref group colums before every group (nColRefEachGrp==ndim)
- Distance to reference: */
- pos = y[1 + nColRefOnce + nColRefEachGrp + g*(nColRefEachGrp+nColPerGrp)][i];
+ pos = y[column][i];
break;
default:
gmx_fatal(FARGS, "Bad error, this error should have been catched before. Ups.\n");
}
}
- /* printf("grp %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
+ /* printf("crd %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
if (bGetMinMax)
{
if (pos < minfound)
{
fprintf(fpcorr, "%g %g\n", k*dt, corr[k]);
}
- fprintf(fpcorr, "&\n");
+ fprintf(fpcorr, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
}
/* esimate integrated correlation time, fitting is too unstable */
printf(" done\n");
if (fpcorr)
{
- ffclose(fpcorr);
+ gmx_ffclose(fpcorr);
}
/* plot IACT along reaction coordinate */
- fp = xvgropen(fn, "Integrated autocorrelation times", "z", "IACT [ps]", opt->oenv);
- fprintf(fp, "@ s0 symbol 1\n@ s0 symbol size 0.5\n@ s0 line linestyle 0\n");
- fprintf(fp, "# WIN tau(gr1) tau(gr2) ...\n");
- for (i = 0; i < nwins; i++)
+ fp = xvgropen(fn, "Integrated autocorrelation times", xlabel, "IACT [ps]", opt->oenv);
+ if (output_env_get_print_xvgr_codes(opt->oenv))
{
- fprintf(fp, "# %3d ", i);
- for (ig = 0; ig < window[i].nPull; ig++)
+ fprintf(fp, "@ s0 symbol 1\n@ s0 symbol size 0.5\n@ s0 line linestyle 0\n");
+ fprintf(fp, "# WIN tau(gr1) tau(gr2) ...\n");
+ for (i = 0; i < nwins; i++)
{
- fprintf(fp, " %11g", window[i].tau[ig]);
+ fprintf(fp, "# %3d ", i);
+ for (ig = 0; ig < window[i].nPull; ig++)
+ {
+ fprintf(fp, " %11g", window[i].tau[ig]);
+ }
+ fprintf(fp, "\n");
}
- fprintf(fp, "\n");
}
for (i = 0; i < nwins; i++)
{
opt->sigSmoothIact);
/* smooth IACT along reaction coordinate and overwrite g=1+2tau */
smoothIact(window, nwins, opt);
- fprintf(fp, "&\n@ s1 symbol 1\n@ s1 symbol size 0.5\n@ s1 line linestyle 0\n");
- fprintf(fp, "@ s1 symbol color 2\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
+ if (output_env_get_print_xvgr_codes(opt->oenv))
+ {
+ fprintf(fp, "@ s1 symbol 1\n@ s1 symbol size 0.5\n@ s1 line linestyle 0\n");
+ fprintf(fp, "@ s1 symbol color 2\n");
+ }
for (i = 0; i < nwins; i++)
{
for (ig = 0; ig < window[i].nPull; ig++)
}
}
}
- ffclose(fp);
+ gmx_ffclose(fp);
printf("Wrote %s\n", fn);
}
nHist++;
fAv += window[i].forceAv[ig];
}
- /* at the same time, rememer closest histogram */
+ /* at the same time, remember closest histogram */
if (dist < distmin)
{
winmin = i;
}
if (opt->verbose)
{
- fp = xvgropen("pmfintegrated.xvg", "PMF from force integration", "z", "PMF [kJ/mol]", opt->oenv);
+ fp = xvgropen("pmfintegrated.xvg", "PMF from force integration", xlabel, "PMF (kJ/mol)", opt->oenv);
for (j = 0; j < opt->bins; ++j)
{
fprintf(fp, "%g %g\n", (j+0.5)*dz+opt->min, pot[j]);
}
- ffclose(fp);
+ gmx_ffclose(fp);
printf("verbose mode: wrote %s with PMF from interated forces\n", "pmfintegrated.xvg");
}
char fmt[1024], fmtign[1024];
int block = 1, sizenow;
- fp = ffopen(opt->fnGroupsel, "r");
+ fp = gmx_ffopen(opt->fnGroupsel, "r");
opt->groupsel = NULL;
snew(tmpbuf, len);
//! Boolean XOR
#define WHAMBOOLXOR(a, b) ( ((!(a)) && (b)) || ((a) && (!(b))))
-/*! Number of elements in fnm (used for command line parsing) */
+//! Number of elements in fnm (used for command line parsing)
#define NFILE asize(fnm)
//! The main g_wham routine
int gmx_wham(int argc, char *argv[])
{
const char *desc[] = {
- "This is an analysis program that implements the Weighted",
+ "[THISMODULE] is an analysis program that implements the Weighted",
"Histogram Analysis Method (WHAM). It is intended to analyze",
"output files generated by umbrella sampling simulations to ",
"compute a potential of mean force (PMF). [PAR] ",
"[TT]*[tt] With option [TT]-ip[tt], the user provides file names of (gzipped) [TT].pdo[tt] files, i.e.",
" the GROMACS 3.3 umbrella output files. If you have some unusual"
" reaction coordinate you may also generate your own [TT].pdo[tt] files and",
- " feed them with the [TT]-ip[tt] option into to [TT]g_wham[tt]. The [TT].pdo[tt] file header",
+ " feed them with the [TT]-ip[tt] option into to [THISMODULE]. The [TT].pdo[tt] file header",
" must be similar to the following:[PAR]",
"[TT]# UMBRELLA 3.0[BR]",
"# Component selection: 0 0 1[BR]",
"If you want the free energy at a different ",
"position to be zero, set [TT]-zprof0[tt] (useful with bootstrapping, see below).[PAR]",
"For cyclic or periodic reaction coordinates (dihedral angle, channel PMF",
- "without osmotic gradient), the option [TT]-cycl[tt] is useful. [TT]g_wham[tt] will make use of the ",
+ "without osmotic gradient), the option [TT]-cycl[tt] is useful.",
+ "[THISMODULE] will make use of the",
"periodicity of the system and generate a periodic PMF. The first and the last bin of the",
"reaction coordinate will assumed be be neighbors.[PAR]",
"Option [TT]-sym[tt] symmetrizes the profile around z=0 before output, ",
"which may be useful for, e.g. membranes.[PAR]",
"AUTOCORRELATIONS[BR]----------------[BR]",
- "With [TT]-ac[tt], [TT]g_wham[tt] estimates the integrated autocorrelation ",
+ "With [TT]-ac[tt], [THISMODULE] estimates the integrated autocorrelation ",
"time (IACT) [GRK]tau[grk] for each umbrella window and weights the respective ",
"window with 1/[1+2*[GRK]tau[grk]/dt]. The IACTs are written ",
"to the file defined with [TT]-oiact[tt]. In verbose mode, all ",
"integration of the ACFs while the ACFs are larger 0.05.",
"If you prefer to compute the IACTs by a more sophisticated (but possibly ",
"less robust) method such as fitting to a double exponential, you can ",
- "compute the IACTs with [TT]g_analyze[tt] and provide them to [TT]g_wham[tt] with the file ",
+ "compute the IACTs with [gmx-analyze] and provide them to [THISMODULE] with the file ",
"[TT]iact-in.dat[tt] (option [TT]-iiact[tt]), which should contain one line per ",
"input file ([TT].pdo[tt] or pullx/f file) and one column per pull group in the respective file.[PAR]",
"ERROR ANALYSIS[BR]--------------[BR]",
opt.stepchange = 100;
opt.stepUpdateContrib = 100;
- if (!parse_common_args(&argc, argv, PCA_BE_NICE,
+ if (!parse_common_args(&argc, argv, 0,
NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &opt.oenv))
{
return 0;
/* write histograms */
histout = xvgropen(opt2fn("-hist", NFILE, fnm), "Umbrella histograms",
- "z", "count", opt.oenv);
+ xlabel, "count", opt.oenv);
for (l = 0; l < opt.bins; ++l)
{
fprintf(histout, "%e\t", (double)(l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
}
fprintf(histout, "\n");
}
- ffclose(histout);
+ gmx_ffclose(histout);
printf("Wrote %s\n", opt2fn("-hist", NFILE, fnm));
if (opt.bHistOnly)
{
}
/* write profile or density of states */
- profout = xvgropen(opt2fn("-o", NFILE, fnm), title, "z", ylabel, opt.oenv);
+ profout = xvgropen(opt2fn("-o", NFILE, fnm), title, xlabel, ylabel, opt.oenv);
for (i = 0; i < opt.bins; ++i)
{
fprintf(profout, "%e\t%e\n", (double)(i+0.5)/opt.bins*(opt.max-opt.min)+opt.min, profile[i]);
}
- ffclose(profout);
+ gmx_ffclose(profout);
printf("Wrote %s\n", opt2fn("-o", NFILE, fnm));
/* Bootstrap Method */
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff