/*! \brief
* enum for energy units
*/
-enum {
- enSel, en_kJ, en_kCal, en_kT, enNr
+enum
+{
+ enSel,
+ en_kJ,
+ en_kCal,
+ en_kT,
+ enNr
};
/*! \brief
* enum for type of input files (pdos, tpr, or pullf)
*/
-enum {
- whamin_unknown, whamin_tpr, whamin_pullxf, whamin_pdo
+enum
+{
+ whamin_unknown,
+ whamin_tpr,
+ whamin_pullxf,
+ whamin_pdo
};
/*! \brief enum for bootstrapping method
* J Chem Theory Comput, 6(12), 3713-3720 (2010)
* ********************************************************************
*/
-enum {
- bsMethod_unknown, bsMethod_BayesianHist, bsMethod_hist,
- bsMethod_traj, bsMethod_trajGauss
+enum
+{
+ bsMethod_unknown,
+ bsMethod_BayesianHist,
+ bsMethod_hist,
+ bsMethod_traj,
+ bsMethod_trajGauss
};
//! Parameters of one pull coodinate
typedef struct
{
- int pull_type; //!< such as constraint, umbrella, ...
- int geometry; //!< such as distance, direction, cylinder
- int ngroup; //!< the number of pull groups involved
- ivec dim; //!< pull dimension with geometry distance
- int ndim; //!< nr of pull_dim != 0
- real k; //!< force constants in tpr file
- real init_dist; //!< reference displacement
- char coord_unit[256]; //!< unit of the displacement
+ int pull_type; //!< such as constraint, umbrella, ...
+ int geometry; //!< such as distance, direction, cylinder
+ int ngroup; //!< the number of pull groups involved
+ ivec dim; //!< pull dimension with geometry distance
+ int ndim; //!< nr of pull_dim != 0
+ real k; //!< force constants in tpr file
+ real init_dist; //!< reference displacement
+ char coord_unit[256]; //!< unit of the displacement
} t_pullcoord;
//! Parameters of the umbrella potentials
*/
/*!\{*/
int npullcrds; //!< nr of umbrella pull coordinates for reading
- t_pullcoord *pcrd; //!< the pull coordinates
+ t_pullcoord* pcrd; //!< the pull coordinates
gmx_bool bPrintCOM; //!< COMs of pull groups writtn in pullx.xvg
gmx_bool bPrintRefValue; //!< Reference value for the coordinate written in pullx.xvg
gmx_bool bPrintComp; //!< Components of pull distance written to pullx.xvg ?
//! Data in the umbrella histograms
typedef struct
{
- int nPull; //!< nr of pull groups in this pdo or pullf/x file
- double **Histo; //!< nPull histograms
- double **cum; //!< nPull cumulative distribution functions
- int nBin; //!< nr of bins. identical to opt->bins
- double *k; //!< force constants for the nPull coords
- double *pos; //!< umbrella positions for the nPull coords
- double *z; //!< z=(-Fi/kT) for the nPull coords. These values are iteratively computed during wham
- int *N; //!< nr of data points in nPull histograms.
- int *Ntot; //!< also nr of data points. N and Ntot only differ if bHistEq==TRUE
+ int nPull; //!< nr of pull groups in this pdo or pullf/x file
+ double** Histo; //!< nPull histograms
+ double** cum; //!< nPull cumulative distribution functions
+ int nBin; //!< nr of bins. identical to opt->bins
+ double* k; //!< force constants for the nPull coords
+ double* pos; //!< umbrella positions for the nPull coords
+ double* z; //!< z=(-Fi/kT) for the nPull coords. These values are iteratively computed during wham
+ int* N; //!< nr of data points in nPull histograms.
+ int* Ntot; //!< also nr of data points. N and Ntot only differ if bHistEq==TRUE
/*! \brief g = 1 + 2*tau[int]/dt where tau is the integrated autocorrelation time.
*
* Compare, e.g. Ferrenberg/Swendsen, PRL 63:1195 (1989),
* Kumar et al, J Comp Chem 13, 1011-1021 (1992), eq. 28
*/
- double *g;
- double *tau; //!< intetrated autocorrelation time (IACT)
- double *tausmooth; //!< smoothed IACT
+ double* g;
+ double* tau; //!< intetrated autocorrelation time (IACT)
+ double* tausmooth; //!< smoothed IACT
- double dt; //!< timestep in the input data. Can be adapted with gmx wham option -dt
+ double dt; //!< timestep in the input data. Can be adapted with gmx wham option -dt
/*! \brief TRUE, if any data point of the histogram is within min and max, otherwise FALSE */
- gmx_bool **bContrib;
- real **ztime; //!< input data z(t) as a function of time. Required to compute ACTs
+ gmx_bool** bContrib;
+ real** ztime; //!< input data z(t) as a function of time. Required to compute ACTs
/*! \brief average force estimated from average displacement, fAv=dzAv*k
*
* Used for integration to guess the potential.
*/
- real *forceAv;
- real *aver; //!< average of histograms
- real *sigma; //!< stddev of histograms
- double *bsWeight; //!< for bootstrapping complete histograms with continuous weights
+ real* forceAv;
+ real* aver; //!< average of histograms
+ real* sigma; //!< stddev of histograms
+ double* bsWeight; //!< for bootstrapping complete histograms with continuous weights
} t_UmbrellaWindow;
//! Selection of pull coordinates to be used in WHAM (one structure for each tpr file)
typedef struct
{
- int n; //!< total nr of pull coords in this tpr file
- int nUse; //!< nr of pull coords used
- gmx_bool *bUse; //!< boolean array of size n. =1 if used, =0 if not
+ int n; //!< total nr of pull coords in this tpr file
+ int nUse; //!< nr of pull coords used
+ gmx_bool* bUse; //!< boolean array of size n. =1 if used, =0 if not
} t_coordselection;
//! Parameters of WHAM
* \name Input stuff
*/
/*!\{*/
- const char *fnTpr, *fnPullf, *fnCoordSel;
- const char *fnPdo, *fnPullx; //!< file names of input
- gmx_bool bTpr, bPullf, bPdo, bPullx; //!< input file types given?
- real tmin, tmax, dt; //!< only read input within tmin and tmax with dt
-
- gmx_bool bInitPotByIntegration; //!< before WHAM, guess potential by force integration. Yields 1.5 to 2 times faster convergence
- int stepUpdateContrib; //!< update contribution table every ... iterations. Accelerates WHAM.
- int nCoordsel; //!< if >0: use only certain group in WHAM, if ==0: use all groups
- t_coordselection *coordsel; //!< for each tpr file: which pull coordinates to use in WHAM?
+ const char *fnTpr, *fnPullf, *fnCoordSel;
+ const char *fnPdo, *fnPullx; //!< file names of input
+ gmx_bool bTpr, bPullf, bPdo, bPullx; //!< input file types given?
+ real tmin, tmax, dt; //!< only read input within tmin and tmax with dt
+
+ gmx_bool bInitPotByIntegration; //!< before WHAM, guess potential by force integration. Yields 1.5 to 2 times faster convergence
+ int stepUpdateContrib; //!< update contribution table every ... iterations. Accelerates WHAM.
+ int nCoordsel; //!< if >0: use only certain group in WHAM, if ==0: use all groups
+ t_coordselection* coordsel; //!< for each tpr file: which pull coordinates to use in WHAM?
/*!\}*/
/*!
* \name Basic WHAM options
*/
/*!\{*/
- int bins; //!< nr of bins, min, max, and dz of profile
+ int bins; //!< nr of bins, min, max, and dz of profile
real min, max, dz;
real Temperature, Tolerance; //!< temperature, converged when probability changes less than Tolerance
gmx_bool bCycl; //!< generate cyclic (periodic) PMF
* \name Output control
*/
/*!\{*/
- gmx_bool bLog; //!< energy output (instead of probability) for profile
- int unit; //!< unit for PMF output kJ/mol or kT or kCal/mol
- gmx_bool bSym; //!< symmetrize PMF around z=0 after WHAM, useful for membranes etc.
+ gmx_bool bLog; //!< energy output (instead of probability) for profile
+ int unit; //!< unit for PMF output kJ/mol or kT or kCal/mol
+ gmx_bool bSym; //!< symmetrize PMF around z=0 after WHAM, useful for membranes etc.
/*! \brief after wham, set prof to zero at this z-position.
* When bootstrapping, set zProf0 to a "stable" reference position.
*/
- real zProf0;
- gmx_bool bProf0Set; //!< setting profile to 0 at zProf0?
+ real zProf0;
+ gmx_bool bProf0Set; //!< setting profile to 0 at zProf0?
- gmx_bool bBoundsOnly, bHistOnly; //!< determine min and max, or write histograms and exit
- gmx_bool bAuto; //!< determine min and max automatically but do not exit
+ gmx_bool bBoundsOnly, bHistOnly; //!< determine min and max, or write histograms and exit
+ gmx_bool bAuto; //!< determine min and max automatically but do not exit
- gmx_bool verbose; //!< more noisy wham mode
- int stepchange; //!< print maximum change in prof after how many interations
- gmx_output_env_t *oenv; //!< xvgr options
+ gmx_bool verbose; //!< more noisy wham mode
+ int stepchange; //!< print maximum change in prof after how many interations
+ gmx_output_env_t* oenv; //!< xvgr options
/*!\}*/
/*!
* \name Autocorrelation stuff
/*!\{*/
gmx_bool bTauIntGiven, bCalcTauInt; //!< IACT given or should be calculated?
real sigSmoothIact; //!< sigma of Gaussian to smooth ACTs
- gmx_bool bAllowReduceIact; //!< Allow to reduce ACTs during smoothing. Otherwise ACT are only increased during smoothing
- real acTrestart; //!< when computing ACT, time between restarting points
+ gmx_bool bAllowReduceIact; //!< Allow to reduce ACTs during smoothing. Otherwise ACT are only increased during smoothing
+ real acTrestart; //!< when computing ACT, time between restarting points
/* \brief Enforce the same weight for each umbella window, that is
* calculate with the same number of data points for
* \name Bootstrapping stuff
*/
/*!\{*/
- int nBootStrap; //!< nr of bootstraps (50 is usually enough)
+ int nBootStrap; //!< nr of bootstraps (50 is usually enough)
/* \brief bootstrap method
*
long the reaction coordinate xi. Avoids gaps along xi. */
int histBootStrapBlockLength;
- int bsSeed; //!< random seed for bootstrapping
+ int bsSeed; //!< random seed for bootstrapping
/* \brief Write cumulative distribution functions (CDFs) of histograms
and write the generated histograms for each bootstrap */
* \name tabulated umbrella potential stuff
*/
/*!\{*/
- gmx_bool bTab;
- double *tabX, *tabY, tabMin, tabMax, tabDz;
- int tabNbins;
+ gmx_bool bTab;
+ double * tabX, *tabY, tabMin, tabMax, tabDz;
+ int tabNbins;
/*!\}*/
- gmx::DefaultRandomEngine rng; //!< gromacs random number generator
- gmx::TabulatedNormalDistribution<> normalDistribution; //!< Uses default: real output, 14-bit table
+ gmx::DefaultRandomEngine rng; //!< gromacs random number generator
+ gmx::TabulatedNormalDistribution<> normalDistribution; //!< Uses default: real output, 14-bit table
} t_UmbrellaOptions;
//! Make an umbrella window (may contain several histograms)
-static t_UmbrellaWindow * initUmbrellaWindows(int nwin)
+static t_UmbrellaWindow* initUmbrellaWindows(int nwin)
{
- t_UmbrellaWindow *win;
+ t_UmbrellaWindow* win;
int i;
snew(win, nwin);
for (i = 0; i < nwin; i++)
{
- win[i].Histo = win[i].cum = nullptr;
- win[i].k = win[i].pos = win[i].z = nullptr;
- win[i].N = win[i].Ntot = nullptr;
- win[i].g = win[i].tau = win[i].tausmooth = nullptr;
- win[i].bContrib = nullptr;
- win[i].ztime = nullptr;
- win[i].forceAv = nullptr;
- win[i].aver = win[i].sigma = nullptr;
- win[i].bsWeight = nullptr;
+ win[i].Histo = win[i].cum = nullptr;
+ win[i].k = win[i].pos = win[i].z = nullptr;
+ win[i].N = win[i].Ntot = nullptr;
+ win[i].g = win[i].tau = win[i].tausmooth = nullptr;
+ win[i].bContrib = nullptr;
+ win[i].ztime = nullptr;
+ win[i].forceAv = nullptr;
+ win[i].aver = win[i].sigma = nullptr;
+ win[i].bsWeight = nullptr;
}
return win;
}
//! Delete an umbrella window (may contain several histograms)
-static void freeUmbrellaWindows(t_UmbrellaWindow *win, int nwin)
+static void freeUmbrellaWindows(t_UmbrellaWindow* win, int nwin)
{
int i, j;
for (i = 0; i < nwin; i++)
/*! \brief
* Read and setup tabulated umbrella potential
*/
-static void setup_tab(const char *fn, t_UmbrellaOptions *opt)
+static void setup_tab(const char* fn, t_UmbrellaOptions* opt)
{
int i, ny, nl;
- double **y;
+ double** y;
printf("Setting up tabulated potential from file %s\n", fn);
nl = read_xvg(fn, &y, &ny);
gmx_fatal(FARGS, "Found %d columns in %s. Expected 2.\n", ny, fn);
}
opt->tabMin = y[0][0];
- opt->tabMax = y[0][nl-1];
- opt->tabDz = (opt->tabMax-opt->tabMin)/(nl-1);
+ opt->tabMax = y[0][nl - 1];
+ opt->tabDz = (opt->tabMax - opt->tabMin) / (nl - 1);
if (opt->tabDz <= 0)
{
- gmx_fatal(FARGS, "The tabulated potential in %s must be provided in \n"
- "ascending z-direction", fn);
+ gmx_fatal(FARGS,
+ "The tabulated potential in %s must be provided in \n"
+ "ascending z-direction",
+ fn);
}
- for (i = 0; i < nl-1; i++)
+ for (i = 0; i < nl - 1; i++)
{
- if (std::abs(y[0][i+1]-y[0][i]-opt->tabDz) > opt->tabDz/1e6)
+ if (std::abs(y[0][i + 1] - y[0][i] - opt->tabDz) > opt->tabDz / 1e6)
{
gmx_fatal(FARGS, "z-values in %s are not equally spaced.\n", fn);
}
opt->tabX[i] = y[0][i];
opt->tabY[i] = y[1][i];
}
- printf("Found equally spaced tabulated potential from %g to %g, spacing %g\n",
- opt->tabMin, opt->tabMax, opt->tabDz);
+ printf("Found equally spaced tabulated potential from %g to %g, spacing %g\n", opt->tabMin,
+ opt->tabMax, opt->tabDz);
}
//! Read the header of an PDO file (position, force const, nr of groups)
-static void read_pdo_header(FILE * file, t_UmbrellaHeader * header, t_UmbrellaOptions *opt)
+static void read_pdo_header(FILE* file, t_UmbrellaHeader* header, t_UmbrellaOptions* opt)
{
char line[2048];
char Buffer0[256], Buffer1[256], Buffer2[256], Buffer3[256], Buffer4[256];
ist >> Buffer0 >> Buffer1 >> Buffer2;
if (std::strcmp(Buffer1, "UMBRELLA") != 0)
{
- gmx_fatal(FARGS, "This does not appear to be a valid pdo file. Found %s, expected %s\n"
+ gmx_fatal(FARGS,
+ "This does not appear to be a valid pdo file. Found %s, expected %s\n"
"(Found in first line: `%s')\n",
Buffer1, "UMBRELLA", line);
}
if (opt->verbose)
{
- printf("\tFound nPull=%d , nSkip=%d, ref=%s\n", header->nPull, header->nSkip,
- header->Reference);
+ printf("\tFound nPull=%d , nSkip=%d, ref=%s\n", header->nPull, header->nSkip, header->Reference);
}
for (i = 0; i < header->nPull; ++i)
if (opt->verbose)
{
- printf("\tpullgroup %d, pullname = %s, UmbPos = %g, UmbConst = %g\n",
- i, header->PullName[i], header->UmbPos[i][0], header->UmbCons[i][0]);
+ printf("\tpullgroup %d, pullname = %s, UmbPos = %g, UmbConst = %g\n", i,
+ header->PullName[i], header->UmbPos[i][0], header->UmbCons[i][0]);
}
}
}
//! smarter fgets
-static char *fgets3(FILE *fp, char ptr[], int *len)
+static char* fgets3(FILE* fp, char ptr[], int* len)
{
- char *p;
+ char* p;
int slen;
- if (fgets(ptr, *len-1, fp) == nullptr)
+ if (fgets(ptr, *len - 1, fp) == nullptr)
{
return nullptr;
}
{
/* This line is longer than len characters, let's increase len! */
*len += STRLEN;
- p += STRLEN;
+ p += STRLEN;
srenew(ptr, *len);
- if (fgets(p-1, STRLEN, fp) == nullptr)
+ if (fgets(p - 1, STRLEN, fp) == nullptr)
{
break;
}
}
slen = std::strlen(ptr);
- if (ptr[slen-1] == '\n')
+ if (ptr[slen - 1] == '\n')
{
- ptr[slen-1] = '\0';
+ ptr[slen - 1] = '\0';
}
return ptr;
* At the moment, this warning is avoided by hiding the format string
* the variable fmtlf.
*/
-static void read_pdo_data(FILE * file, t_UmbrellaHeader * header,
- int fileno, t_UmbrellaWindow * win,
- t_UmbrellaOptions *opt,
- gmx_bool bGetMinMax, real *mintmp, real *maxtmp)
+static void read_pdo_data(FILE* file,
+ t_UmbrellaHeader* header,
+ int fileno,
+ t_UmbrellaWindow* win,
+ t_UmbrellaOptions* opt,
+ gmx_bool bGetMinMax,
+ real* mintmp,
+ real* maxtmp)
{
- int i, inttemp, bins, count, ntot;
- real minval, maxval, minfound = 1e20, maxfound = -1e20;
- double temp, time, time0 = 0, dt;
- char *ptr = nullptr;
- t_UmbrellaWindow * window = nullptr;
- gmx_bool timeok, dt_ok = true;
- char *tmpbuf = nullptr, fmt[256], fmtign[256], fmtlf[5] = "%lf";
- int len = STRLEN, dstep = 1;
- const int blocklen = 4096;
- int *lennow = nullptr;
+ int i, inttemp, bins, count, ntot;
+ real minval, maxval, minfound = 1e20, maxfound = -1e20;
+ double temp, time, time0 = 0, dt;
+ char* ptr = nullptr;
+ t_UmbrellaWindow* window = nullptr;
+ gmx_bool timeok, dt_ok = true;
+ char * tmpbuf = nullptr, fmt[256], fmtign[256], fmtlf[5] = "%lf";
+ int len = STRLEN, dstep = 1;
+ const int blocklen = 4096;
+ int* lennow = nullptr;
if (!bGetMinMax)
{
- bins = opt->bins;
- minval = opt->min;
- maxval = opt->max;
+ bins = opt->bins;
+ minval = opt->min;
+ maxval = opt->max;
- window = win+fileno;
+ window = win + fileno;
/* Need to alocate memory and set up structure */
window->nPull = header->nPull;
window->nBin = bins;
{
minfound = 1e20;
maxfound = -1e20;
- minval = maxval = bins = 0; /* Get rid of warnings */
+ minval = maxval = bins = 0; /* Get rid of warnings */
}
count = 0;
snew(tmpbuf, len);
- while ( (ptr = fgets3(file, tmpbuf, &len)) != nullptr)
+ while ((ptr = fgets3(file, tmpbuf, &len)) != nullptr)
{
trim(ptr);
sscanf(ptr, fmtlf, &time); /* printf("Time %f\n",time); */
/* Round time to fs */
- time = 1.0/1000*( gmx::roundToInt64(time*1000) );
+ time = 1.0 / 1000 * (gmx::roundToInt64(time * 1000));
/* get time step of pdo file */
if (count == 0)
}
else if (count == 1)
{
- dt = time-time0;
+ dt = time - time0;
if (opt->dt > 0.0)
{
- dstep = gmx::roundToInt(opt->dt/dt);
+ dstep = gmx::roundToInt(opt->dt / dt);
if (dstep == 0)
{
dstep = 1;
}
if (!bGetMinMax)
{
- window->dt = dt*dstep;
+ window->dt = dt * dstep;
}
}
count++;
- dt_ok = ((count-1)%dstep == 0);
+ dt_ok = ((count - 1) % dstep == 0);
timeok = (dt_ok && time >= opt->tmin && time <= opt->tmax);
/* if (opt->verbose)
printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
for (i = 0; i < header->nPull; ++i)
{
std::strcpy(fmt, fmtign);
- std::strcat(fmt, "%lf"); /* Creating a format stings such as "%*s...%*s%lf" */
- std::strcat(fmtign, "%*s"); /* ignoring one more entry in the next loop */
+ std::strcat(fmt, "%lf"); /* Creating a format stings such as "%*s...%*s%lf" */
+ std::strcat(fmtign, "%*s"); /* ignoring one more entry in the next loop */
if (sscanf(ptr, fmt, &temp))
{
temp += header->UmbPos[i][0];
}
temp -= minval;
- temp /= (maxval-minval);
+ temp /= (maxval - minval);
temp *= bins;
- temp = std::floor(temp);
+ temp = std::floor(temp);
- inttemp = static_cast<int> (temp);
+ inttemp = static_cast<int>(temp);
if (opt->bCycl)
{
if (inttemp < 0)
* by this routine (not recommended). Since we now support autocorrelations, it is better to set
* an appropriate autocorrelation times instead of using this function.
*/
-static void enforceEqualWeights(t_UmbrellaWindow * window, int nWindows)
+static void enforceEqualWeights(t_UmbrellaWindow* window, int nWindows)
{
int i, k, j, NEnforced;
double ratio;
NEnforced = window[0].Ntot[0];
printf("\nFound -hist-eq. Enforcing equal weights for all histograms, \ni.e. doing a "
- "non-weighted histogram analysis method. Ndata = %d\n", NEnforced);
+ "non-weighted histogram analysis method. Ndata = %d\n",
+ NEnforced);
/* enforce all histograms to have the same weight as the very first histogram */
for (j = 0; j < nWindows; ++j)
{
for (k = 0; k < window[j].nPull; ++k)
{
- ratio = 1.0*NEnforced/window[j].Ntot[k];
+ ratio = 1.0 * NEnforced / window[j].Ntot[k];
for (i = 0; i < window[0].nBin; ++i)
{
window[j].Histo[k][i] *= ratio;
}
- window[j].N[k] = gmx::roundToInt(ratio*window[j].N[k]);
+ window[j].N[k] = gmx::roundToInt(ratio * window[j].N[k]);
}
}
}
/*! \brief Simple linear interpolation between two given tabulated points
*/
-static double tabulated_pot(double dist, t_UmbrellaOptions *opt)
+static double tabulated_pot(double dist, t_UmbrellaOptions* opt)
{
int jl, ju;
double pl, pu, dz, dp;
- jl = static_cast<int> (std::floor((dist-opt->tabMin)/opt->tabDz));
- ju = jl+1;
+ jl = static_cast<int>(std::floor((dist - opt->tabMin) / opt->tabDz));
+ ju = jl + 1;
if (jl < 0 || ju >= opt->tabNbins)
{
- gmx_fatal(FARGS, "Distance %f out of bounds of tabulated potential (jl=%d, ju=%d).\n"
- "Provide an extended table.", dist, jl, ju);
+ gmx_fatal(FARGS,
+ "Distance %f out of bounds of tabulated potential (jl=%d, ju=%d).\n"
+ "Provide an extended table.",
+ dist, jl, ju);
}
pl = opt->tabY[jl];
pu = opt->tabY[ju];
- dz = dist-opt->tabX[jl];
- dp = (pu-pl)*dz/opt->tabDz;
- return pl+dp;
+ dz = dist - opt->tabX[jl];
+ dp = (pu - pl) * dz / opt->tabDz;
+ return pl + dp;
}
* After rapid convergence (using only substiantal contributions), we always switch to
* full precision.
*/
-static void setup_acc_wham(const double *profile, t_UmbrellaWindow * window, int nWindows,
- t_UmbrellaOptions *opt)
+static void setup_acc_wham(const double* profile, t_UmbrellaWindow* window, int nWindows, t_UmbrellaOptions* opt)
{
- int i, j, k, nGrptot = 0, nContrib = 0, nTot = 0;
- double U, min = opt->min, dz = opt->dz, temp, ztot_half, distance, ztot, contrib1, contrib2;
- gmx_bool bAnyContrib;
- static int bFirst = 1;
+ int i, j, k, nGrptot = 0, nContrib = 0, nTot = 0;
+ double U, min = opt->min, dz = opt->dz, temp, ztot_half, distance, ztot, contrib1, contrib2;
+ gmx_bool bAnyContrib;
+ static int bFirst = 1;
static double wham_contrib_lim;
if (bFirst)
{
nGrptot += window[i].nPull;
}
- wham_contrib_lim = opt->Tolerance/nGrptot;
+ wham_contrib_lim = opt->Tolerance / nGrptot;
}
- ztot = opt->max-opt->min;
- ztot_half = ztot/2;
+ ztot = opt->max - opt->min;
+ ztot_half = ztot / 2;
for (i = 0; i < nWindows; ++i)
{
bAnyContrib = FALSE;
for (k = 0; k < opt->bins; ++k)
{
- temp = (1.0*k+0.5)*dz+min;
- distance = temp - window[i].pos[j]; /* distance to umbrella center */
+ temp = (1.0 * k + 0.5) * dz + min;
+ distance = temp - window[i].pos[j]; /* distance to umbrella center */
if (opt->bCycl)
- { /* in cyclic wham: */
- if (distance > ztot_half) /* |distance| < ztot_half */
+ { /* in cyclic wham: */
+ if (distance > ztot_half) /* |distance| < ztot_half */
{
distance -= ztot;
}
if (!opt->bTab)
{
- U = 0.5*window[i].k[j]*gmx::square(distance); /* harmonic potential assumed. */
+ U = 0.5 * window[i].k[j] * gmx::square(distance); /* harmonic potential assumed. */
}
else
{
- U = tabulated_pot(distance, opt); /* Use tabulated potential */
+ U = tabulated_pot(distance, opt); /* Use tabulated potential */
}
- contrib1 = profile[k]*std::exp(-U/(BOLTZ*opt->Temperature));
- contrib2 = window[i].N[j]*std::exp(-U/(BOLTZ*opt->Temperature) + window[i].z[j]);
+ contrib1 = profile[k] * std::exp(-U / (BOLTZ * opt->Temperature));
+ contrib2 = window[i].N[j] * std::exp(-U / (BOLTZ * opt->Temperature) + window[i].z[j]);
window[i].bContrib[j][k] = (contrib1 > wham_contrib_lim || contrib2 > wham_contrib_lim);
bAnyContrib = bAnyContrib || window[i].bContrib[j][k];
if (window[i].bContrib[j][k])
if (bFirst)
{
printf("Initialized rapid wham stuff (contrib tolerance %g)\n"
- "Evaluating only %d of %d expressions.\n\n", wham_contrib_lim, nContrib, nTot);
+ "Evaluating only %d of %d expressions.\n\n",
+ wham_contrib_lim, nContrib, nTot);
}
if (opt->verbose)
{
- printf("Updated rapid wham stuff. (evaluating only %d of %d contributions)\n",
- nContrib, nTot);
+ printf("Updated rapid wham stuff. (evaluating only %d of %d contributions)\n", nContrib, nTot);
}
bFirst = 0;
}
//! Compute the PMF (one of the two main WHAM routines)
-static void calc_profile(double *profile, t_UmbrellaWindow * window, int nWindows,
- t_UmbrellaOptions *opt, gmx_bool bExact)
+static void calc_profile(double* profile, t_UmbrellaWindow* window, int nWindows, t_UmbrellaOptions* opt, gmx_bool bExact)
{
double ztot_half, ztot, min = opt->min, dz = opt->dz;
- ztot = opt->max-opt->min;
- ztot_half = ztot/2;
+ ztot = opt->max - opt->min;
+ ztot_half = ztot / 2;
#pragma omp parallel
{
int nthreads = gmx_omp_get_max_threads();
int thread_id = gmx_omp_get_thread_num();
int i;
- int i0 = thread_id*opt->bins/nthreads;
- int i1 = std::min(opt->bins, ((thread_id+1)*opt->bins)/nthreads);
+ int i0 = thread_id * opt->bins / nthreads;
+ int i1 = std::min(opt->bins, ((thread_id + 1) * opt->bins) / nthreads);
for (i = i0; i < i1; ++i)
{
{
for (k = 0; k < window[j].nPull; ++k)
{
- invg = 1.0/window[j].g[k] * window[j].bsWeight[k];
- temp = (1.0*i+0.5)*dz+min;
- num += invg*window[j].Histo[k][i];
+ invg = 1.0 / window[j].g[k] * window[j].bsWeight[k];
+ temp = (1.0 * i + 0.5) * dz + min;
+ num += invg * window[j].Histo[k][i];
if (!(bExact || window[j].bContrib[k][i]))
{
continue;
}
- distance = temp - window[j].pos[k]; /* distance to umbrella center */
+ distance = temp - window[j].pos[k]; /* distance to umbrella center */
if (opt->bCycl)
- { /* in cyclic wham: */
- if (distance > ztot_half) /* |distance| < ztot_half */
+ { /* in cyclic wham: */
+ if (distance > ztot_half) /* |distance| < ztot_half */
{
distance -= ztot;
}
if (!opt->bTab)
{
- U = 0.5*window[j].k[k]*gmx::square(distance); /* harmonic potential assumed. */
+ U = 0.5 * window[j].k[k] * gmx::square(distance); /* harmonic potential assumed. */
}
else
{
- U = tabulated_pot(distance, opt); /* Use tabulated potential */
+ U = tabulated_pot(distance, opt); /* Use tabulated potential */
}
- denom += invg*window[j].N[k]*std::exp(-U/(BOLTZ*opt->Temperature) + window[j].z[k]);
+ denom += invg * window[j].N[k]
+ * std::exp(-U / (BOLTZ * opt->Temperature) + window[j].z[k]);
}
}
- profile[i] = num/denom;
+ profile[i] = num / denom;
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
}
//! Compute the free energy offsets z (one of the two main WHAM routines)
-static double calc_z(const double * profile, t_UmbrellaWindow * window, int nWindows,
- t_UmbrellaOptions *opt, gmx_bool bExact)
+static double calc_z(const double* profile, t_UmbrellaWindow* window, int nWindows, t_UmbrellaOptions* opt, gmx_bool bExact)
{
- double min = opt->min, dz = opt->dz, ztot_half, ztot;
+ double min = opt->min, dz = opt->dz, ztot_half, ztot;
double maxglob = -1e20;
- ztot = opt->max-opt->min;
- ztot_half = ztot/2;
+ ztot = opt->max - opt->min;
+ ztot_half = ztot / 2;
#pragma omp parallel
{
int nthreads = gmx_omp_get_max_threads();
int thread_id = gmx_omp_get_thread_num();
int i;
- int i0 = thread_id*nWindows/nthreads;
- int i1 = std::min(nWindows, ((thread_id+1)*nWindows)/nthreads);
- double maxloc = -1e20;
+ int i0 = thread_id * nWindows / nthreads;
+ int i1 = std::min(nWindows, ((thread_id + 1) * nWindows) / nthreads);
+ double maxloc = -1e20;
for (i = i0; i < i1; ++i)
{
- double total = 0, temp, distance, U = 0;
+ double total = 0, temp, distance, U = 0;
int j, k;
for (j = 0; j < window[i].nPull; ++j)
{
continue;
}
- temp = (1.0*k+0.5)*dz+min;
- distance = temp - window[i].pos[j]; /* distance to umbrella center */
+ temp = (1.0 * k + 0.5) * dz + min;
+ distance = temp - window[i].pos[j]; /* distance to umbrella center */
if (opt->bCycl)
- { /* in cyclic wham: */
- if (distance > ztot_half) /* |distance| < ztot_half */
+ { /* in cyclic wham: */
+ if (distance > ztot_half) /* |distance| < ztot_half */
{
distance -= ztot;
}
if (!opt->bTab)
{
- U = 0.5*window[i].k[j]*gmx::square(distance); /* harmonic potential assumed. */
+ U = 0.5 * window[i].k[j] * gmx::square(distance); /* harmonic potential assumed. */
}
else
{
- U = tabulated_pot(distance, opt); /* Use tabulated potential */
+ U = tabulated_pot(distance, opt); /* Use tabulated potential */
}
- total += profile[k]*std::exp(-U/(BOLTZ*opt->Temperature));
+ total += profile[k] * std::exp(-U / (BOLTZ * opt->Temperature));
}
/* Avoid floating point exception if window is far outside min and max */
if (total != 0.0)
}
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
return maxglob;
}
//! Make PMF symmetric around 0 (useful e.g. for membranes)
-static void symmetrizeProfile(double* profile, t_UmbrellaOptions *opt)
+static void symmetrizeProfile(double* profile, t_UmbrellaOptions* opt)
{
int i, j, bins = opt->bins;
double *prof2, min = opt->min, max = opt->max, dz = opt->dz, zsym, deltaz, profsym;
if (min > 0. || max < 0.)
{
- gmx_fatal(FARGS, "Cannot symmetrize profile around z=0 with min=%f and max=%f\n",
- opt->min, opt->max);
+ gmx_fatal(FARGS, "Cannot symmetrize profile around z=0 with min=%f and max=%f\n", opt->min,
+ opt->max);
}
snew(prof2, bins);
for (i = 0; i < bins; i++)
{
- z = min+(i+0.5)*dz;
+ z = min + (i + 0.5) * dz;
zsym = -z;
/* bin left of zsym */
- j = gmx::roundToInt((zsym-min)/dz)-1;
- if (j >= 0 && (j+1) < bins)
+ j = gmx::roundToInt((zsym - min) / dz) - 1;
+ if (j >= 0 && (j + 1) < bins)
{
/* interpolate profile linearly between bins j and j+1 */
- z1 = min+(j+0.5)*dz;
- deltaz = zsym-z1;
- profsym = profile[j] + (profile[j+1]-profile[j])/dz*deltaz;
+ z1 = min + (j + 0.5) * dz;
+ deltaz = zsym - z1;
+ profsym = profile[j] + (profile[j + 1] - profile[j]) / dz * deltaz;
/* average between left and right */
- prof2[i] = 0.5*(profsym+profile[i]);
+ prof2[i] = 0.5 * (profsym + profile[i]);
}
else
{
}
}
- std::memcpy(profile, prof2, bins*sizeof(double));
+ std::memcpy(profile, prof2, bins * sizeof(double));
sfree(prof2);
}
//! Set energy unit (kJ/mol,kT,kCal/mol) and set it to zero at opt->zProf0
-static void prof_normalization_and_unit(double * profile, t_UmbrellaOptions *opt)
+static void prof_normalization_and_unit(double* profile, t_UmbrellaOptions* opt)
{
int i, bins, imin;
double unit_factor = 1., diff;
- bins = opt->bins;
+ bins = opt->bins;
/* Not log? Nothing to do! */
if (!opt->bLog)
}
else if (opt->unit == en_kJ)
{
- unit_factor = BOLTZ*opt->Temperature;
+ unit_factor = BOLTZ * opt->Temperature;
}
else if (opt->unit == en_kCal)
{
- unit_factor = (BOLTZ/CAL2JOULE)*opt->Temperature;
+ unit_factor = (BOLTZ / CAL2JOULE) * opt->Temperature;
}
else
{
{
if (profile[i] > 0.0)
{
- profile[i] = -std::log(profile[i])*unit_factor;
+ profile[i] = -std::log(profile[i]) * unit_factor;
}
}
{
/* Get bin with shortest distance to opt->zProf0
(-0.5 from bin position and +0.5 from rounding cancel) */
- imin = static_cast<int>((opt->zProf0-opt->min)/opt->dz);
+ imin = static_cast<int>((opt->zProf0 - opt->min) / opt->dz);
if (imin < 0)
{
imin = 0;
}
else if (imin >= bins)
{
- imin = bins-1;
+ imin = bins - 1;
}
diff = profile[imin];
}
}
//! Make an array of random integers (used for bootstrapping)
-static void getRandomIntArray(int nPull, int blockLength, int* randomArray, gmx::DefaultRandomEngine * rng)
+static void getRandomIntArray(int nPull, int blockLength, int* randomArray, gmx::DefaultRandomEngine* rng)
{
- gmx::UniformIntDistribution<int> dist(0, blockLength-1);
+ gmx::UniformIntDistribution<int> dist(0, blockLength - 1);
int ipull, blockBase, nr, ipullRandom;
for (ipull = 0; ipull < nPull; ipull++)
{
- blockBase = (ipull/blockLength)*blockLength;
+ blockBase = (ipull / blockLength) * blockLength;
do
{ /* make sure nothing bad happens in the last block */
nr = dist(*rng); // [0,blockLength-1]
ipullRandom = blockBase + nr;
- }
- while (ipullRandom >= nPull);
+ } while (ipullRandom >= nPull);
if (ipullRandom < 0 || ipullRandom >= nPull)
{
- gmx_fatal(FARGS, "Ups, random iWin = %d, nPull = %d, nr = %d, "
+ gmx_fatal(FARGS,
+ "Ups, random iWin = %d, nPull = %d, nr = %d, "
"blockLength = %d, blockBase = %d\n",
ipullRandom, nPull, nr, blockLength, blockBase);
}
* This is used when bootstapping new trajectories and thereby create new histogtrams,
* but it is not required if we bootstrap complete histograms.
*/
-static void copy_pullgrp_to_synthwindow(t_UmbrellaWindow *synthWindow,
- t_UmbrellaWindow *thisWindow, int pullid)
+static void copy_pullgrp_to_synthwindow(t_UmbrellaWindow* synthWindow, t_UmbrellaWindow* thisWindow, int pullid)
{
- synthWindow->N [0] = thisWindow->N [pullid];
- synthWindow->Histo [0] = thisWindow->Histo [pullid];
- synthWindow->pos [0] = thisWindow->pos [pullid];
- synthWindow->z [0] = thisWindow->z [pullid];
- synthWindow->k [0] = thisWindow->k [pullid];
- synthWindow->bContrib[0] = thisWindow->bContrib [pullid];
- synthWindow->g [0] = thisWindow->g [pullid];
- synthWindow->bsWeight[0] = thisWindow->bsWeight [pullid];
+ synthWindow->N[0] = thisWindow->N[pullid];
+ synthWindow->Histo[0] = thisWindow->Histo[pullid];
+ synthWindow->pos[0] = thisWindow->pos[pullid];
+ synthWindow->z[0] = thisWindow->z[pullid];
+ synthWindow->k[0] = thisWindow->k[pullid];
+ synthWindow->bContrib[0] = thisWindow->bContrib[pullid];
+ synthWindow->g[0] = thisWindow->g[pullid];
+ synthWindow->bsWeight[0] = thisWindow->bsWeight[pullid];
}
/*! \brief Calculate cumulative distribution function of of all histograms.
* which are distributed according to the histograms. Required to generate
* the "synthetic" histograms for the Bootstrap method
*/
-static void calc_cumulatives(t_UmbrellaWindow *window, int nWindows,
- t_UmbrellaOptions *opt, const char *fnhist, const char *xlabel)
+static void calc_cumulatives(t_UmbrellaWindow* window,
+ int nWindows,
+ t_UmbrellaOptions* opt,
+ const char* fnhist,
+ const char* xlabel)
{
int i, j, k, nbin;
double last;
std::string fn;
- FILE *fp = nullptr;
+ FILE* fp = nullptr;
if (opt->bs_verbose)
{
snew(window[i].cum, window[i].nPull);
for (j = 0; j < window[i].nPull; j++)
{
- snew(window[i].cum[j], nbin+1);
+ snew(window[i].cum[j], nbin + 1);
window[i].cum[j][0] = 0.;
for (k = 1; k <= nbin; k++)
{
- window[i].cum[j][k] = window[i].cum[j][k-1]+window[i].Histo[j][k-1];
+ window[i].cum[j][k] = window[i].cum[j][k - 1] + window[i].Histo[j][k - 1];
}
/* normalize CDFs. Ensure cum[nbin]==1 */
{
for (k = 0; k <= nbin; k++)
{
- fprintf(fp, "%g\t", opt->min+k*opt->dz);
+ fprintf(fp, "%g\t", opt->min + k * opt->dz);
for (i = 0; i < nWindows; i++)
{
for (j = 0; j < window[i].nPull; j++)
*
* This is used to generate a random sequence distributed according to a histogram
*/
-static void searchCumulative(const double xx[], int n, double x, int *j)
+static void searchCumulative(const double xx[], int n, double x, int* j)
{
int ju, jm, jl;
jl = -1;
ju = n;
- while (ju-jl > 1)
+ while (ju - jl > 1)
{
- jm = (ju+jl) >> 1;
+ jm = (ju + jl) >> 1;
if (x >= xx[jm])
{
jl = jm;
{
*j = 0;
}
- else if (x == xx[n-1])
+ else if (x == xx[n - 1])
{
- *j = n-2;
+ *j = n - 2;
}
else
{
}
//! Bootstrap new trajectories and thereby generate new (bootstrapped) histograms
-static void create_synthetic_histo(t_UmbrellaWindow *synthWindow, t_UmbrellaWindow *thisWindow,
- int pullid, t_UmbrellaOptions *opt)
+static void create_synthetic_histo(t_UmbrellaWindow* synthWindow,
+ t_UmbrellaWindow* thisWindow,
+ int pullid,
+ t_UmbrellaOptions* opt)
{
int N, i, nbins, r_index, ibin;
double r, tausteps = 0.0, a, ap, dt, x, invsqrt2, g, y, sig = 0., z, mu = 0.;
/* tau = autocorrelation time */
if (opt->tauBootStrap > 0.0)
{
- tausteps = opt->tauBootStrap/dt;
+ tausteps = opt->tauBootStrap / dt;
}
else if (opt->bTauIntGiven || opt->bCalcTauInt)
{
/* calc tausteps from g=1+2tausteps */
g = thisWindow->g[pullid];
- tausteps = (g-1)/2;
+ tausteps = (g - 1) / 2;
}
else
{
gmx_fatal(FARGS, "%s", errstr);
}
- synthWindow->N [0] = N;
- synthWindow->pos [0] = thisWindow->pos[pullid];
- synthWindow->z [0] = thisWindow->z[pullid];
- synthWindow->k [0] = thisWindow->k[pullid];
+ synthWindow->N[0] = N;
+ synthWindow->pos[0] = thisWindow->pos[pullid];
+ synthWindow->z[0] = thisWindow->z[pullid];
+ synthWindow->k[0] = thisWindow->k[pullid];
synthWindow->bContrib[0] = thisWindow->bContrib[pullid];
- synthWindow->g [0] = thisWindow->g [pullid];
+ synthWindow->g[0] = thisWindow->g[pullid];
synthWindow->bsWeight[0] = thisWindow->bsWeight[pullid];
for (i = 0; i < nbins; i++)
if (opt->bsMethod == bsMethod_trajGauss)
{
- sig = thisWindow->sigma [pullid];
- mu = thisWindow->aver [pullid];
+ sig = thisWindow->sigma[pullid];
+ mu = thisWindow->aver[pullid];
}
/* Genrate autocorrelated Gaussian random variable with autocorrelation time tau
Note: The ACT of the flat distribution and of the generated histogram is not
100% exactly tau, but near tau (my test was 3.8 instead of 4).
*/
- a = std::exp(-1.0/tausteps);
- ap = std::sqrt(1.0-a*a);
- invsqrt2 = 1.0/std::sqrt(2.0);
+ a = std::exp(-1.0 / tausteps);
+ ap = std::sqrt(1.0 - a * a);
+ invsqrt2 = 1.0 / std::sqrt(2.0);
/* init random sequence */
x = opt->normalDistribution(opt->rng);
for (i = 0; i < N; i++)
{
y = opt->normalDistribution(opt->rng);
- x = a*x+ap*y;
+ x = a * x + ap * y;
/* get flat distribution in [0,1] using cumulative distribution function of Gauusian
Note: CDF(Gaussian) = 0.5*{1+erf[x/sqrt(2)]}
*/
- r = 0.5*(1+std::erf(x*invsqrt2));
- searchCumulative(thisWindow->cum[pullid], nbins+1, r, &r_index);
+ r = 0.5 * (1 + std::erf(x * invsqrt2));
+ searchCumulative(thisWindow->cum[pullid], nbins + 1, r, &r_index);
synthWindow->Histo[0][r_index] += 1.;
}
}
while (i < N)
{
y = opt->normalDistribution(opt->rng);
- x = a*x+ap*y;
- z = x*sig+mu;
- ibin = static_cast<int> (std::floor((z-opt->min)/opt->dz));
+ x = a * x + ap * y;
+ z = x * sig + mu;
+ ibin = static_cast<int>(std::floor((z - opt->min) / opt->dz));
if (opt->bCycl)
{
if (ibin < 0)
{
- while ( (ibin += nbins) < 0)
- {
- ;
- }
+ while ((ibin += nbins) < 0) {}
}
else if (ibin >= nbins)
{
- while ( (ibin -= nbins) >= nbins)
- {
- ;
- }
+ while ((ibin -= nbins) >= nbins) {}
}
}
* If bs_index>=0, a number is added to the output file name to allow the ouput of all
* sets of bootstrapped histograms.
*/
-static void print_histograms(const char *fnhist, t_UmbrellaWindow * window, int nWindows,
- int bs_index, t_UmbrellaOptions *opt, const char *xlabel)
+static void print_histograms(const char* fnhist,
+ t_UmbrellaWindow* window,
+ int nWindows,
+ int bs_index,
+ t_UmbrellaOptions* opt,
+ const char* xlabel)
{
std::string fn, title;
- FILE *fp;
+ FILE* fp;
int bins, l, i, j;
if (bs_index >= 0)
/* Write histograms */
for (l = 0; l < bins; ++l)
{
- fprintf(fp, "%e\t", (l+0.5)*opt->dz+opt->min);
+ fprintf(fp, "%e\t", (l + 0.5) * opt->dz + opt->min);
for (i = 0; i < nWindows; ++i)
{
for (j = 0; j < window[i].nPull; ++j)
}
//! Make random weights for histograms for the Bayesian bootstrap of complete histograms)
-static void setRandomBsWeights(t_UmbrellaWindow *synthwin, int nAllPull, t_UmbrellaOptions *opt)
+static void setRandomBsWeights(t_UmbrellaWindow* synthwin, int nAllPull, t_UmbrellaOptions* opt)
{
- int i;
- double *r;
+ int i;
+ double* r;
gmx::UniformRealDistribution<real> dist(0, nAllPull);
snew(r, nAllPull);
/* generate ordered random numbers between 0 and nAllPull */
- for (i = 0; i < nAllPull-1; i++)
+ for (i = 0; i < nAllPull - 1; i++)
{
r[i] = dist(opt->rng);
}
- std::sort(r, r+nAllPull-1);
- r[nAllPull-1] = 1.0*nAllPull;
+ std::sort(r, r + nAllPull - 1);
+ r[nAllPull - 1] = 1.0 * nAllPull;
synthwin[0].bsWeight[0] = r[0];
for (i = 1; i < nAllPull; i++)
{
- synthwin[i].bsWeight[0] = r[i]-r[i-1];
+ synthwin[i].bsWeight[0] = r[i] - r[i - 1];
}
/* avoid to have zero weight by adding a tiny value */
}
//! The main bootstrapping routine
-static void do_bootstrapping(const char *fnres, const char* fnprof, const char *fnhist,
- const char *xlabel, char* ylabel, double *profile,
- t_UmbrellaWindow * window, int nWindows, t_UmbrellaOptions *opt)
+static void do_bootstrapping(const char* fnres,
+ const char* fnprof,
+ const char* fnhist,
+ const char* xlabel,
+ char* ylabel,
+ double* profile,
+ t_UmbrellaWindow* window,
+ int nWindows,
+ t_UmbrellaOptions* opt)
{
- t_UmbrellaWindow * synthWindow;
- double *bsProfile, *bsProfiles_av, *bsProfiles_av2, maxchange = 1e20, tmp, stddev;
- int i, j, *randomArray = nullptr, winid, pullid, ib;
- int iAllPull, nAllPull, *allPull_winId, *allPull_pullId;
- FILE *fp;
- gmx_bool bExact = FALSE;
+ t_UmbrellaWindow* synthWindow;
+ double * bsProfile, *bsProfiles_av, *bsProfiles_av2, maxchange = 1e20, tmp, stddev;
+ int i, j, *randomArray = nullptr, winid, pullid, ib;
+ int iAllPull, nAllPull, *allPull_winId, *allPull_pullId;
+ FILE* fp;
+ gmx_bool bExact = FALSE;
/* init random generator */
if (opt->bsSeed == 0)
}
opt->rng.seed(opt->bsSeed);
- snew(bsProfile, opt->bins);
+ snew(bsProfile, opt->bins);
snew(bsProfiles_av, opt->bins);
snew(bsProfiles_av2, opt->bins);
/* just copy all histogams into synthWindow array */
for (i = 0; i < nAllPull; i++)
{
- winid = allPull_winId [i];
+ winid = allPull_winId[i];
pullid = allPull_pullId[i];
- copy_pullgrp_to_synthwindow(synthWindow+i, window+winid, pullid);
+ copy_pullgrp_to_synthwindow(synthWindow + i, window + winid, pullid);
}
break;
case bsMethod_traj:
- case bsMethod_trajGauss:
- calc_cumulatives(window, nWindows, opt, fnhist, xlabel);
- break;
- default:
- gmx_fatal(FARGS, "Unknown bootstrap method. That should not have happened.\n");
+ case bsMethod_trajGauss: calc_cumulatives(window, nWindows, opt, fnhist, xlabel); break;
+ default: gmx_fatal(FARGS, "Unknown bootstrap method. That should not have happened.\n");
}
/* do bootstrapping */
{
printf(" *******************************************\n"
" ******** Start bootstrap nr %d ************\n"
- " *******************************************\n", ib+1);
+ " *******************************************\n",
+ ib + 1);
switch (opt->bsMethod)
{
getRandomIntArray(nAllPull, opt->histBootStrapBlockLength, randomArray, &opt->rng);
for (i = 0; i < nAllPull; i++)
{
- winid = allPull_winId [randomArray[i]];
+ winid = allPull_winId[randomArray[i]];
pullid = allPull_pullId[randomArray[i]];
- copy_pullgrp_to_synthwindow(synthWindow+i, window+winid, pullid);
+ copy_pullgrp_to_synthwindow(synthWindow + i, window + winid, pullid);
}
break;
case bsMethod_BayesianHist:
{
winid = allPull_winId[i];
pullid = allPull_pullId[i];
- create_synthetic_histo(synthWindow+i, window+winid, pullid, opt);
+ create_synthetic_histo(synthWindow + i, window + winid, pullid, opt);
}
break;
}
i = 0;
bExact = FALSE;
maxchange = 1e20;
- std::memcpy(bsProfile, profile, opt->bins*sizeof(double)); /* use profile as guess */
+ std::memcpy(bsProfile, profile, opt->bins * sizeof(double)); /* use profile as guess */
do
{
- if ( (i%opt->stepUpdateContrib) == 0)
+ if ((i % opt->stepUpdateContrib) == 0)
{
setup_acc_wham(bsProfile, synthWindow, nAllPull, opt);
}
{
bExact = TRUE;
}
- if (((i%opt->stepchange) == 0 || i == 1) && i != 0)
+ if (((i % opt->stepchange) == 0 || i == 1) && i != 0)
{
printf("\t%4d) Maximum change %e\n", i, maxchange);
}
calc_profile(bsProfile, synthWindow, nAllPull, opt, bExact);
i++;
- }
- while ( (maxchange = calc_z(bsProfile, synthWindow, nAllPull, opt, bExact)) > opt->Tolerance || !bExact);
+ } while ((maxchange = calc_z(bsProfile, synthWindow, nAllPull, opt, bExact)) > opt->Tolerance
+ || !bExact);
printf("\tConverged in %d iterations. Final maximum change %g\n", i, maxchange);
if (opt->bLog)
/* save stuff to get average and stddev */
for (i = 0; i < opt->bins; i++)
{
- tmp = bsProfile[i];
- bsProfiles_av[i] += tmp;
- bsProfiles_av2[i] += tmp*tmp;
- fprintf(fp, "%e\t%e\n", (i+0.5)*opt->dz+opt->min, tmp);
+ tmp = bsProfile[i];
+ bsProfiles_av[i] += tmp;
+ bsProfiles_av2[i] += tmp * tmp;
+ fprintf(fp, "%e\t%e\n", (i + 0.5) * opt->dz + opt->min, tmp);
}
fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
}
}
for (i = 0; i < opt->bins; i++)
{
- bsProfiles_av [i] /= opt->nBootStrap;
+ bsProfiles_av[i] /= opt->nBootStrap;
bsProfiles_av2[i] /= opt->nBootStrap;
- tmp = bsProfiles_av2[i]-gmx::square(bsProfiles_av[i]);
- stddev = (tmp >= 0.) ? std::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);
+ tmp = bsProfiles_av2[i] - gmx::square(bsProfiles_av[i]);
+ stddev = (tmp >= 0.) ? std::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);
}
xvgrclose(fp);
printf("Wrote boot strap result to %s\n", fnres);
}
//! Return type of input file based on file extension (xvg, pdo, or tpr)
-static int whaminFileType(char *fn)
+static int whaminFileType(char* fn)
{
int len;
len = std::strlen(fn);
- if (std::strcmp(fn+len-3, "tpr") == 0)
+ if (std::strcmp(fn + len - 3, "tpr") == 0)
{
return whamin_tpr;
}
- else if (std::strcmp(fn+len-3, "xvg") == 0 || std::strcmp(fn+len-6, "xvg.gz") == 0)
+ else if (std::strcmp(fn + len - 3, "xvg") == 0 || std::strcmp(fn + len - 6, "xvg.gz") == 0)
{
return whamin_pullxf;
}
- else if (std::strcmp(fn+len-3, "pdo") == 0 || std::strcmp(fn+len-6, "pdo.gz") == 0)
+ else if (std::strcmp(fn + len - 3, "pdo") == 0 || std::strcmp(fn + len - 6, "pdo.gz") == 0)
{
return whamin_pdo;
}
}
//! Read the files names in pdo-files.dat, pullf/x-files.dat, tpr-files.dat
-static void read_wham_in(const char *fn, char ***filenamesRet, int *nfilesRet,
- t_UmbrellaOptions *opt)
+static void read_wham_in(const char* fn, char*** filenamesRet, int* nfilesRet, t_UmbrellaOptions* opt)
{
- char **filename = nullptr, tmp[WHAM_MAXFILELEN+2];
+ char **filename = nullptr, tmp[WHAM_MAXFILELEN + 2];
int nread, sizenow, i, block = 1;
- FILE *fp;
+ FILE* fp;
fp = gmx_ffopen(fn, "r");
nread = 0;
{
sizenow += block;
srenew(filename, sizenow);
- for (i = sizenow-block; i < sizenow; i++)
+ for (i = sizenow - block; i < sizenow; i++)
{
snew(filename[i], WHAM_MAXFILELEN);
}
}
/* remove newline if there is one */
- if (tmp[std::strlen(tmp)-1] == '\n')
+ if (tmp[std::strlen(tmp) - 1] == '\n')
{
- tmp[std::strlen(tmp)-1] = '\0';
+ tmp[std::strlen(tmp) - 1] = '\0';
}
std::strcpy(filename[nread], tmp);
if (opt->verbose)
}
//! Open a file or a pipe to a gzipped file
-static FILE *open_pdo_pipe(const char *fn, t_UmbrellaOptions *opt, gmx_bool *bPipeOpen)
+static FILE* open_pdo_pipe(const char* fn, t_UmbrellaOptions* opt, gmx_bool* bPipeOpen)
{
char Buffer[2048], gunzip[1024], *Path = nullptr;
- FILE *pipe = nullptr;
+ FILE* pipe = nullptr;
static gmx_bool bFirst = true;
/* gzipped pdo file? */
- if ((std::strcmp(fn+std::strlen(fn)-3, ".gz") == 0))
+ if ((std::strcmp(fn + std::strlen(fn) - 3, ".gz") == 0))
{
/* search gunzip executable */
if (!(Path = getenv("GMX_PATH_GZIP")))
}
else
{
- gmx_fatal(FARGS, "Cannot find executable gunzip in /bin or /usr/bin.\n"
+ gmx_fatal(FARGS,
+ "Cannot find executable gunzip in /bin or /usr/bin.\n"
"You may want to define the path to gunzip "
"with the environment variable GMX_PATH_GZIP.");
}
sprintf(gunzip, "%s/gunzip", Path);
if (!gmx_fexist(gunzip))
{
- gmx_fatal(FARGS, "Cannot find executable %s. Please define the path to gunzip"
- " in the environmental varialbe GMX_PATH_GZIP.", gunzip);
+ gmx_fatal(FARGS,
+ "Cannot find executable %s. Please define the path to gunzip"
+ " in the environmental varialbe GMX_PATH_GZIP.",
+ gunzip);
}
}
if (bFirst)
}
//! Close file or pipe
-static void pdo_close_file(FILE *fp)
+static void pdo_close_file(FILE* fp)
{
#if HAVE_PIPES
pclose(fp);
}
//! Reading all pdo files
-static void read_pdo_files(char **fn, int nfiles, t_UmbrellaHeader* header,
- t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
+static void read_pdo_files(char** fn, int nfiles, t_UmbrellaHeader* header, t_UmbrellaWindow* window, t_UmbrellaOptions* opt)
{
- FILE *file;
+ FILE* file;
real mintmp, maxtmp, done = 0.;
int i;
gmx_bool bPipeOpen;
file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
/*fgets(Buffer0,999,file);
fprintf(stderr,"First line '%s'\n",Buffer0); */
- done = 100.0*(i+1)/nfiles;
- fprintf(stdout, "\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
+ done = 100.0 * (i + 1) / nfiles;
+ fprintf(stdout, "\rOpening %s ... [%2.0f%%]", fn[i], done);
+ fflush(stdout);
if (opt->verbose)
{
printf("\n");
if (opt->bBoundsOnly)
{
printf("Found option -boundsonly, now exiting.\n");
- exit (0);
+ exit(0);
}
}
/* store stepsize in profile */
- opt->dz = (opt->max-opt->min)/opt->bins;
+ opt->dz = (opt->max - opt->min) / opt->bins;
/* Having min and max, we read in all files */
/* Loop over all files */
for (i = 0; i < nfiles; ++i)
{
- done = 100.0*(i+1)/nfiles;
- fprintf(stdout, "\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
+ done = 100.0 * (i + 1) / nfiles;
+ fprintf(stdout, "\rOpening %s ... [%2.0f%%]", fn[i], done);
+ fflush(stdout);
if (opt->verbose)
{
printf("\n");
read_pdo_header(file, header, opt);
/* load data into window */
read_pdo_data(file, header, i, window, opt, FALSE, nullptr, nullptr);
- if ((window+i)->Ntot[0] == 0)
+ if ((window + i)->Ntot[0] == 0)
{
fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fn[i]);
}
#define int2YN(a) (((a) == 0) ? ("N") : ("Y"))
//! Read pull groups from a tpr file (including position, force const, geometry, number of groups)
-static void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions *opt, t_coordselection *coordsel)
+static void read_tpr_header(const char* fn, t_UmbrellaHeader* header, t_UmbrellaOptions* opt, t_coordselection* coordsel)
{
- t_inputrec irInstance;
- t_inputrec *ir = &irInstance;
- t_state state;
- static int first = 1;
+ t_inputrec irInstance;
+ t_inputrec* ir = &irInstance;
+ t_state state;
+ static int first = 1;
/* printf("Reading %s \n",fn); */
read_tpx_state(fn, ir, &state, nullptr);
snew(header->pcrd, header->npullcrds);
for (int i = 0; i < ir->pull->ncoord; i++)
{
- header->pcrd[i].pull_type = ir->pull->coord[i].eType;
- header->pcrd[i].geometry = ir->pull->coord[i].eGeom;
- header->pcrd[i].ngroup = ir->pull->coord[i].ngroup;
+ header->pcrd[i].pull_type = ir->pull->coord[i].eType;
+ header->pcrd[i].geometry = ir->pull->coord[i].eGeom;
+ header->pcrd[i].ngroup = ir->pull->coord[i].ngroup;
/* Angle type coordinates are handled fully in degrees in gmx wham.
* The pull "distance" appears with a power of -2 in the force constant,
* so we need to multiply with the internal units (radians for angle)
* to user units (degrees for an angle) with the same power.
*/
- header->pcrd[i].k = ir->pull->coord[i].k/gmx::square(pull_conversion_factor_internal2userinput(&ir->pull->coord[i]));
- header->pcrd[i].init_dist = ir->pull->coord[i].init;
+ header->pcrd[i].k =
+ ir->pull->coord[i].k
+ / gmx::square(pull_conversion_factor_internal2userinput(&ir->pull->coord[i]));
+ header->pcrd[i].init_dist = ir->pull->coord[i].init;
copy_ivec(ir->pull->coord[i].dim, header->pcrd[i].dim);
- header->pcrd[i].ndim = header->pcrd[i].dim[XX] + header->pcrd[i].dim[YY] + header->pcrd[i].dim[ZZ];
+ header->pcrd[i].ndim =
+ header->pcrd[i].dim[XX] + header->pcrd[i].dim[YY] + header->pcrd[i].dim[ZZ];
- std::strcpy(header->pcrd[i].coord_unit,
- pull_coordinate_units(&ir->pull->coord[i]));
+ std::strcpy(header->pcrd[i].coord_unit, pull_coordinate_units(&ir->pull->coord[i]));
if (ir->efep != efepNO && ir->pull->coord[i].k != ir->pull->coord[i].kB)
{
- gmx_fatal(FARGS, "Seems like you did free-energy perturbation, and you perturbed the force constant."
+ gmx_fatal(FARGS,
+ "Seems like you did free-energy perturbation, and you perturbed the force "
+ "constant."
" This is not supported.\n");
}
if (coordsel && (coordsel->n != ir->pull->ncoord))
{
- gmx_fatal(FARGS, "Found %d pull coordinates in %s, but %d columns in the respective line\n"
- "coordinate selection file (option -is)\n", ir->pull->ncoord, fn, coordsel->n);
+ gmx_fatal(FARGS,
+ "Found %d pull coordinates in %s, but %d columns in the respective line\n"
+ "coordinate selection file (option -is)\n",
+ ir->pull->ncoord, fn, coordsel->n);
}
}
{
if (header->pcrd[i].pull_type != epullUMBRELLA)
{
- gmx_fatal(FARGS, "%s: Pull coordinate %d is of type \"%s\", expected \"umbrella\". Only umbrella coodinates can enter WHAM.\n"
- "If you have umrella and non-umbrella coordinates, you can select the umbrella coordinates with gmx wham -is\n",
- fn, i+1, epull_names[header->pcrd[i].pull_type]);
+ gmx_fatal(FARGS,
+ "%s: Pull coordinate %d is of type \"%s\", expected \"umbrella\". Only "
+ "umbrella coodinates can enter WHAM.\n"
+ "If you have umrella and non-umbrella coordinates, you can select the "
+ "umbrella coordinates with gmx wham -is\n",
+ fn, i + 1, epull_names[header->pcrd[i].pull_type]);
}
if (!geometryIsSet)
{
}
if (geom != header->pcrd[i].geometry)
{
- gmx_fatal(FARGS, "%s: Your pull coordinates have different pull geometry (coordinate 1: %s, coordinate %d: %s)\n"
- "If you want to use only some pull coordinates in WHAM, please select them with option gmx wham -is\n",
- fn, epullg_names[geom], i+1, epullg_names[header->pcrd[i].geometry]);
- }
- if (thedim[XX] != header->pcrd[i].dim[XX] || thedim[YY] != header->pcrd[i].dim[YY] || thedim[ZZ] != header->pcrd[i].dim[ZZ])
- {
- gmx_fatal(FARGS, "%s: Your pull coordinates have different pull dimensions (coordinate 1: %s %s %s, coordinate %d: %s %s %s)\n"
- "If you want to use only some pull coordinates in WHAM, please select them with option gmx wham -is\n",
- fn, int2YN(thedim[XX]), int2YN(thedim[YY]), int2YN(thedim[ZZ]), i+1,
- int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]), int2YN(header->pcrd[i].dim[ZZ]));
+ gmx_fatal(FARGS,
+ "%s: Your pull coordinates have different pull geometry (coordinate 1: "
+ "%s, coordinate %d: %s)\n"
+ "If you want to use only some pull coordinates in WHAM, please select "
+ "them with option gmx wham -is\n",
+ fn, epullg_names[geom], i + 1, epullg_names[header->pcrd[i].geometry]);
+ }
+ if (thedim[XX] != header->pcrd[i].dim[XX] || thedim[YY] != header->pcrd[i].dim[YY]
+ || thedim[ZZ] != header->pcrd[i].dim[ZZ])
+ {
+ gmx_fatal(FARGS,
+ "%s: Your pull coordinates have different pull dimensions (coordinate 1: "
+ "%s %s %s, coordinate %d: %s %s %s)\n"
+ "If you want to use only some pull coordinates in WHAM, please select "
+ "them with option gmx wham -is\n",
+ fn, int2YN(thedim[XX]), int2YN(thedim[YY]), int2YN(thedim[ZZ]), i + 1,
+ int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]),
+ int2YN(header->pcrd[i].dim[ZZ]));
}
if (header->pcrd[i].geometry == epullgCYL)
{
if (header->pcrd[i].dim[XX] || header->pcrd[i].dim[YY] || (!header->pcrd[i].dim[ZZ]))
{
- gmx_fatal(FARGS, "With pull geometry 'cylinder', expected pulling in Z direction only.\n"
- "However, found dimensions [%s %s %s]\n",
- int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]),
- int2YN(header->pcrd[i].dim[ZZ]));
+ gmx_fatal(
+ FARGS,
+ "With pull geometry 'cylinder', expected pulling in Z direction only.\n"
+ "However, found dimensions [%s %s %s]\n",
+ int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]),
+ int2YN(header->pcrd[i].dim[ZZ]));
}
}
if (header->pcrd[i].k <= 0.0)
{
- gmx_fatal(FARGS, "%s: Pull coordinate %d has force constant of of %g.\n"
+ gmx_fatal(FARGS,
+ "%s: Pull coordinate %d has force constant of of %g.\n"
"That doesn't seem to be an Umbrella tpr.\n",
- fn, i+1, header->pcrd[i].k);
+ fn, i + 1, header->pcrd[i].k);
}
}
}
maxlen = (lentmp > maxlen) ? lentmp : maxlen;
}
char fmt[STRLEN];
- sprintf(fmt, "\tGeometry %%-%ds k = %%-8g position = %%-8g dimensions [%%s %%s %%s] (%%d dimensions). Used: %%s\n",
- maxlen+1);
+ sprintf(fmt,
+ "\tGeometry %%-%ds k = %%-8g position = %%-8g dimensions [%%s %%s %%s] (%%d "
+ "dimensions). Used: %%s\n",
+ maxlen + 1);
for (int i = 0; i < ir->pull->ncoord; i++)
{
bool use = (coordsel == nullptr || coordsel->bUse[i]);
- printf(fmt,
- epullg_names[header->pcrd[i].geometry], header->pcrd[i].k, header->pcrd[i].init_dist,
- int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]), int2YN(header->pcrd[i].dim[ZZ]),
- header->pcrd[i].ndim, use ? "Yes" : "No");
- printf("\tPull group coordinates of %d groups expected in pullx files.\n", ir->pull->bPrintCOM ? header->pcrd[i].ngroup : 0);
+ printf(fmt, epullg_names[header->pcrd[i].geometry], header->pcrd[i].k, header->pcrd[i].init_dist,
+ int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]),
+ int2YN(header->pcrd[i].dim[ZZ]), header->pcrd[i].ndim, use ? "Yes" : "No");
+ printf("\tPull group coordinates of %d groups expected in pullx files.\n",
+ ir->pull->bPrintCOM ? header->pcrd[i].ngroup : 0);
}
printf("\tReference value of the coordinate%s expected in pullx files.\n",
header->bPrintRefValue ? "" : " not");
}
//! Read pullx.xvg or pullf.xvg
-static void read_pull_xf(const char *fn, t_UmbrellaHeader * header,
- t_UmbrellaWindow * window,
- t_UmbrellaOptions *opt,
- gmx_bool bGetMinMax, real *mintmp, real *maxtmp,
- t_coordselection *coordsel)
+static void read_pull_xf(const char* fn,
+ t_UmbrellaHeader* header,
+ t_UmbrellaWindow* window,
+ t_UmbrellaOptions* opt,
+ gmx_bool bGetMinMax,
+ real* mintmp,
+ real* maxtmp,
+ t_coordselection* coordsel)
{
- double **y = nullptr, pos = 0., t, force, time0 = 0., dt;
+ double ** y = nullptr, pos = 0., t, force, time0 = 0., dt;
int ny, nt, bins, ibin, i, g, gUsed, dstep = 1;
int nColExpect, ntot, column;
real min, max, minfound = 1e20, maxfound = -1e20;
gmx_bool dt_ok, timeok;
- const char *quantity;
+ const char* quantity;
const int blocklen = 4096;
- int *lennow = nullptr;
+ int* lennow = nullptr;
static gmx_bool bFirst = TRUE;
/*
if (header->bPrintComp && opt->bPullx)
{
- gmx_fatal(FARGS, "gmx wham cannot read pullx files if the components of the coordinate was written\n"
- "(mdp option pull-print-components). Provide the pull force files instead (with option -if).\n");
+ gmx_fatal(
+ FARGS,
+ "gmx wham cannot read pullx files if the components of the coordinate was written\n"
+ "(mdp option pull-print-components). Provide the pull force files instead (with "
+ "option -if).\n");
}
int *nColThisCrd, *nColCOMCrd, *nColRefCrd;
snew(nColThisCrd, header->npullcrds);
- snew(nColCOMCrd, header->npullcrds);
- snew(nColRefCrd, header->npullcrds);
+ snew(nColCOMCrd, header->npullcrds);
+ snew(nColRefCrd, header->npullcrds);
if (!opt->bPullx)
{
/* pullx reading. Note explanation above. */
for (g = 0; g < header->npullcrds; g++)
{
- nColRefCrd[g] = (header->bPrintRefValue ? 1 : 0);
- nColCOMCrd[g] = (header->bPrintCOM ? header->pcrd[g].ndim*header->pcrd[g].ngroup : 0);
- nColThisCrd[g] = 1 + nColCOMCrd[g] + nColRefCrd[g];
+ nColRefCrd[g] = (header->bPrintRefValue ? 1 : 0);
+ nColCOMCrd[g] = (header->bPrintCOM ? header->pcrd[g].ndim * header->pcrd[g].ngroup : 0);
+ nColThisCrd[g] = 1 + nColCOMCrd[g] + nColRefCrd[g];
}
}
nt = read_xvg(fn, &y, &ny);
/* Check consistency */
- quantity = opt->bPullx ? "position" : "force";
+ quantity = opt->bPullx ? "position" : "force";
if (nt < 1)
{
gmx_fatal(FARGS, "Empty pull %s file %s\n", quantity, fn);
printf("\nReading pull %s file %s, expecting %d columns:\n", quantity, fn, nColExpect);
for (i = 0; i < header->npullcrds; i++)
{
- printf("\tColumns for pull coordinate %d\n", i+1);
+ printf("\tColumns for pull coordinate %d\n", i + 1);
printf("\t\treaction coordinate: %d\n"
"\t\tcenter-of-mass of groups: %d\n"
"\t\treference position column: %s\n",
}
if (nColExpect != ny)
{
- gmx_fatal(FARGS, "Expected %d columns (including time column) in %s, but found %d."
- " Maybe you confused options -if and -ix ?", nColExpect, fn, ny);
+ gmx_fatal(FARGS,
+ "Expected %d columns (including time column) in %s, but found %d."
+ " Maybe you confused options -if and -ix ?",
+ nColExpect, fn, ny);
}
if (!bGetMinMax)
max = opt->max;
if (nt > 1)
{
- window->dt = y[0][1]-y[0][0];
+ window->dt = y[0][1] - y[0][0];
}
else if (opt->nBootStrap && opt->tauBootStrap != 0.0)
{
/* Use only groups selected with option -is file */
if (header->npullcrds != coordsel->n)
{
- gmx_fatal(FARGS, "tpr file contains %d pull groups, but expected %d from group selection file\n",
+ gmx_fatal(FARGS,
+ "tpr file contains %d pull groups, but expected %d from group selection "
+ "file\n",
header->npullcrds, coordsel->n);
}
window->nPull = coordsel->nUse;
}
window->nBin = bins;
- snew(window->Histo, window->nPull);
- snew(window->z, window->nPull);
- snew(window->k, window->nPull);
- snew(window->pos, window->nPull);
- snew(window->N, window->nPull);
- snew(window->Ntot, window->nPull);
- snew(window->g, window->nPull);
+ snew(window->Histo, window->nPull);
+ snew(window->z, window->nPull);
+ snew(window->k, window->nPull);
+ snew(window->pos, window->nPull);
+ snew(window->N, window->nPull);
+ snew(window->Ntot, window->nPull);
+ snew(window->g, window->nPull);
snew(window->bsWeight, window->nPull);
window->bContrib = nullptr;
for (g = 0; g < window->nPull; ++g)
{
- window->z [g] = 1;
- window->bsWeight [g] = 1.;
- window->N [g] = 0;
- window->Ntot [g] = 0;
- window->g [g] = 1.;
+ window->z[g] = 1;
+ window->bsWeight[g] = 1.;
+ window->N[g] = 0;
+ window->Ntot[g] = 0;
+ window->g[g] = 1.;
snew(window->Histo[g], bins);
if (opt->bCalcTauInt)
{
continue;
}
- window->k [gUsed] = header->pcrd[g].k;
+ window->k[gUsed] = header->pcrd[g].k;
window->pos[gUsed] = header->pcrd[g].init_dist;
gUsed++;
}
}
else
- { /* only determine min and max */
+ { /* only determine min and max */
minfound = 1e20;
maxfound = -1e20;
- min = max = bins = 0; /* Get rid of warnings */
+ min = max = bins = 0; /* Get rid of warnings */
}
for (i = 0; i < nt; i++)
{
/* Do you want that time frame? */
- t = 1.0/1000*( gmx::roundToInt64((y[0][i]*1000))); /* round time to fs */
+ t = 1.0 / 1000 * (gmx::roundToInt64((y[0][i] * 1000))); /* round time to fs */
/* get time step of pdo file and get dstep from opt->dt */
if (i == 0)
}
else if (i == 1)
{
- dt = t-time0;
+ dt = t - time0;
if (opt->dt > 0.0)
{
- dstep = gmx::roundToInt(opt->dt/dt);
+ dstep = gmx::roundToInt(opt->dt / dt);
if (dstep == 0)
{
dstep = 1;
}
if (!bGetMinMax)
{
- window->dt = dt*dstep;
+ window->dt = dt * dstep;
}
}
- dt_ok = (i%dstep == 0);
+ dt_ok = (i % dstep == 0);
timeok = (dt_ok && t >= opt->tmin && t <= opt->tmax);
/*if (opt->verbose)
printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
* if coordsel != NULL:
* only groups with coordsel.bUse[g]==TRUE are stored. gUsed is not always equal g
*/
- gUsed = -1;
+ gUsed = -1;
for (g = 0; g < header->npullcrds; ++g)
{
/* was this group selected for application in WHAM? */
if (opt->bPullf)
{
/* y has 1 time column y[0] and one column per force y[1],...,y[nCrds] */
- force = y[g+1][i];
- pos = -force/header->pcrd[g].k + header->pcrd[g].init_dist;
+ force = y[g + 1][i];
+ pos = -force / header->pcrd[g].k + header->pcrd[g].init_dist;
}
else
{
{
column += nColThisCrd[j];
}
- pos = y[column][i];
+ pos = y[column][i];
}
/* printf("crd %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
{
if (gUsed >= window->nPull)
{
- gmx_fatal(FARGS, "gUsed too large (%d, nPull=%d). This error should have been caught before.\n",
+ gmx_fatal(FARGS,
+ "gUsed too large (%d, nPull=%d). This error should have been "
+ "caught before.\n",
gUsed, window->nPull);
}
window->ztime[gUsed][ntot] = pos;
}
- ibin = static_cast<int> (std::floor((pos-min)/(max-min)*bins));
+ ibin = static_cast<int>(std::floor((pos - min) / (max - min) * bins));
if (opt->bCycl)
{
if (ibin < 0)
{
- while ( (ibin += bins) < 0)
- {
- ;
- }
+ while ((ibin += bins) < 0) {}
}
else if (ibin >= bins)
{
- while ( (ibin -= bins) >= bins)
- {
- ;
- }
+ while ((ibin -= bins) >= bins) {}
}
}
if (ibin >= 0 && ibin < bins)
}
//! read pullf-files.dat or pullx-files.dat and tpr-files.dat
-static void read_tpr_pullxf_files(char **fnTprs, char **fnPull, int nfiles,
- t_UmbrellaHeader* header,
- t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
+static void read_tpr_pullxf_files(char** fnTprs,
+ char** fnPull,
+ int nfiles,
+ t_UmbrellaHeader* header,
+ t_UmbrellaWindow* window,
+ t_UmbrellaOptions* opt)
{
int i;
real mintmp, maxtmp;
read_tpr_header(fnTprs[i], header, opt, (opt->nCoordsel > 0) ? &opt->coordsel[i] : nullptr);
if (whaminFileType(fnPull[i]) != whamin_pullxf)
{
- gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
+ gmx_fatal(FARGS,
+ "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
}
read_pull_xf(fnPull[i], header, nullptr, opt, TRUE, &mintmp, &maxtmp,
(opt->nCoordsel > 0) ? &opt->coordsel[i] : nullptr);
if (opt->bBoundsOnly)
{
printf("Found option -boundsonly, now exiting.\n");
- exit (0);
+ exit(0);
}
}
/* store stepsize in profile */
- opt->dz = (opt->max-opt->min)/opt->bins;
+ opt->dz = (opt->max - opt->min) / opt->bins;
bool foundData = false;
for (i = 0; i < nfiles; i++)
read_tpr_header(fnTprs[i], header, opt, (opt->nCoordsel > 0) ? &opt->coordsel[i] : nullptr);
if (whaminFileType(fnPull[i]) != whamin_pullxf)
{
- gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
+ gmx_fatal(FARGS,
+ "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
}
- read_pull_xf(fnPull[i], header, window+i, opt, FALSE, nullptr, nullptr,
+ read_pull_xf(fnPull[i], header, window + i, opt, FALSE, nullptr, nullptr,
(opt->nCoordsel > 0) ? &opt->coordsel[i] : nullptr);
if (window[i].Ntot[0] == 0)
{
}
if (!foundData)
{
- gmx_fatal(FARGS, "No data points were found in pullf/pullx files. Maybe you need to specify the -b option?\n");
+ gmx_fatal(FARGS,
+ "No data points were found in pullf/pullx files. Maybe you need to specify the "
+ "-b option?\n");
}
for (i = 0; i < nfiles; i++)
* Note: Here we consider tau[int] := int_0^inf ACF(t) as the integrated autocorrelation times.
* The factor `g := 1 + 2*tau[int]` subsequently enters the uncertainty.
*/
-static void readIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins, const char* fn)
+static void readIntegratedAutocorrelationTimes(t_UmbrellaWindow* window, int nwins, const char* fn)
{
int nlines, ny, i, ig;
- double **iact;
+ double** iact;
printf("Readging Integrated autocorrelation times from %s ...\n", fn);
nlines = read_xvg(fn, &iact, &ny);
{
if (window[i].nPull != ny)
{
- gmx_fatal(FARGS, "You are providing autocorrelation times with option -iiact and the\n"
+ gmx_fatal(FARGS,
+ "You are providing autocorrelation times with option -iiact and the\n"
"number of pull groups is different in different simulations. That is not\n"
"supported yet. Sorry.\n");
}
for (ig = 0; ig < window[i].nPull; ig++)
{
/* compare Kumar et al, J Comp Chem 13, 1011-1021 (1992) */
- window[i].g[ig] = 1+2*iact[ig][i]/window[i].dt;
+ window[i].g[ig] = 1 + 2 * iact[ig][i] / window[i].dt;
if (iact[ig][i] <= 0.0)
{
* If opt->bAllowReduceIact==FALSE, the ACTs are never reduced, only increased
* by the smoothing
*/
-static void smoothIact(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt)
+static void smoothIact(t_UmbrellaWindow* window, int nwins, t_UmbrellaOptions* opt)
{
int i, ig, j, jg;
double pos, dpos2, siglim, siglim2, gaufact, invtwosig2, w, weight, tausm;
/* only evaluate within +- 3sigma of the Gausian */
- siglim = 3.0*opt->sigSmoothIact;
+ siglim = 3.0 * opt->sigSmoothIact;
siglim2 = gmx::square(siglim);
/* pre-factor of Gaussian */
- gaufact = 1.0/(std::sqrt(2*M_PI)*opt->sigSmoothIact);
- invtwosig2 = 0.5/gmx::square(opt->sigSmoothIact);
+ gaufact = 1.0 / (std::sqrt(2 * M_PI) * opt->sigSmoothIact);
+ invtwosig2 = 0.5 / gmx::square(opt->sigSmoothIact);
for (i = 0; i < nwins; i++)
{
{
for (jg = 0; jg < window[j].nPull; jg++)
{
- dpos2 = gmx::square(window[j].pos[jg]-pos);
+ dpos2 = gmx::square(window[j].pos[jg] - pos);
if (dpos2 < siglim2)
{
- w = gaufact*std::exp(-dpos2*invtwosig2);
+ w = gaufact * std::exp(-dpos2 * invtwosig2);
weight += w;
- tausm += w*window[j].tau[jg];
+ tausm += w * window[j].tau[jg];
/*printf("Weight %g dpos2=%g pos=%g gaufact=%g invtwosig2=%g\n",
w,dpos2,pos,gaufact,invtwosig2); */
}
{
window[i].tausmooth[ig] = window[i].tau[ig];
}
- window[i].g[ig] = 1+2*tausm/window[i].dt;
+ window[i].g[ig] = 1 + 2 * tausm / window[i].dt;
}
}
}
/*! \brief Try to compute the autocorrelation time for each umbrealla window
*/
-static void calcIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins,
- t_UmbrellaOptions *opt, const char *fn, const char *xlabel)
+static void calcIntegratedAutocorrelationTimes(t_UmbrellaWindow* window,
+ int nwins,
+ t_UmbrellaOptions* opt,
+ const char* fn,
+ const char* xlabel)
{
int i, ig, ncorr, ntot, j, k, *count, restart;
real *corr, c0, dt, tmp;
printf("\n");
for (i = 0; i < nwins; i++)
{
- fprintf(stdout, "\rEstimating integrated autocorrelation times ... [%2.0f%%] ...", 100.*(i+1)/nwins);
+ fprintf(stdout, "\rEstimating integrated autocorrelation times ... [%2.0f%%] ...",
+ 100. * (i + 1) / nwins);
fflush(stdout);
ntot = window[i].Ntot[0];
/* using half the maximum time as length of autocorrelation function */
- ncorr = ntot/2;
+ ncorr = ntot / 2;
if (ntot < 10)
{
- gmx_fatal(FARGS, "Tryig to estimtate autocorrelation time from only %d"
- " points. Provide more pull data!", ntot);
+ gmx_fatal(FARGS,
+ "Tryig to estimtate autocorrelation time from only %d"
+ " points. Provide more pull data!",
+ ntot);
}
snew(corr, ncorr);
/* snew(corrSq,ncorr); */
snew(count, ncorr);
dt = window[i].dt;
snew(window[i].tau, window[i].nPull);
- restart = gmx::roundToInt(opt->acTrestart/dt);
+ restart = gmx::roundToInt(opt->acTrestart / dt);
if (restart == 0)
{
restart = 1;
{
if (ntot != window[i].Ntot[ig])
{
- gmx_fatal(FARGS, "Encountered different nr of frames in different pull groups.\n"
- "That should not happen. (%d and %d)\n", ntot, window[i].Ntot[ig]);
+ gmx_fatal(FARGS,
+ "Encountered different nr of frames in different pull groups.\n"
+ "That should not happen. (%d and %d)\n",
+ ntot, window[i].Ntot[ig]);
}
ztime = window[i].ztime[ig];
}
for (j = 0; (j < ntot); j += restart)
{
- for (k = 0; (k < ncorr) && (j+k < ntot); k++)
+ for (k = 0; (k < ncorr) && (j + k < ntot); k++)
{
- tmp = (ztime[j]-av)*(ztime[j+k]-av);
- corr [k] += tmp;
+ tmp = (ztime[j] - av) * (ztime[j + k] - av);
+ corr[k] += tmp;
/* corrSq[k] += tmp*tmp; */
count[k]++;
}
for (k = 0; (k < ncorr); k++)
{
/* count probably = (ncorr-k+(restart-1))/restart; */
- corr[k] = corr[k]/count[k];
+ corr[k] = corr[k] / count[k];
/* variance of autocorrelation function */
/* corrSq[k]=corrSq[k]/count[k]; */
}
/* normalize such that corr[0] == 0 */
- c0 = 1./corr[0];
+ c0 = 1. / corr[0];
for (k = 0; (k < ncorr); k++)
{
corr[k] *= c0;
{
for (k = 0; (k < ncorr); k++)
{
- fprintf(fpcorr, "%g %g\n", k*dt, corr[k]);
+ fprintf(fpcorr, "%g %g\n", k * dt, corr[k]);
}
fprintf(fpcorr, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
}
/* esimate integrated correlation time, fitting is too unstable */
- tausteps = 0.5*corr[0];
+ tausteps = 0.5 * corr[0];
/* consider corr below WHAM_AC_ZERO_LIMIT as noise */
for (j = 1; (j < ncorr) && (corr[j] > WHAM_AC_ZERO_LIMIT); j++)
{
/* g = 1+2*tau, see. Ferrenberg/Swendsen, PRL 63:1195 (1989) or
Kumar et al, eq. 28 ff. */
- window[i].tau[ig] = tausteps*dt;
- window[i].g[ig] = 1+2*tausteps;
+ window[i].tau[ig] = tausteps * dt;
+ window[i].g[ig] = 1 + 2 * tausteps;
/* printf("win %d, group %d, estimated correlation time = %g ps\n",i,ig,window[i].tau[ig]); */
} /* ig loop */
sfree(corr);
}
if (opt->sigSmoothIact > 0.0)
{
- printf("Smoothing autocorrelation times along reaction coordinate with Gaussian of sig = %g\n",
+ printf("Smoothing autocorrelation times along reaction coordinate with Gaussian of sig = "
+ "%g\n",
opt->sigSmoothIact);
/* smooth IACT along reaction coordinate and overwrite g=1+2tau */
smoothIact(window, nwins, opt);
/*! \brief
* compute average and sigma of each umbrella histogram
*/
-static void averageSigma(t_UmbrellaWindow *window, int nwins)
+static void averageSigma(t_UmbrellaWindow* window, int nwins)
{
int i, ig, ntot, k;
real av, sum2, sig, diff, *ztime, nSamplesIndep;
av /= ntot;
for (k = 0, sum2 = 0.; k < ntot; k++)
{
- diff = ztime[k]-av;
- sum2 += diff*diff;
+ diff = ztime[k] - av;
+ sum2 += diff * diff;
}
- sig = std::sqrt(sum2/ntot);
+ sig = std::sqrt(sum2 / ntot);
window[i].aver[ig] = av;
/* Note: This estimate for sigma is biased from the limited sampling.
*/
if (window[i].tau)
{
- nSamplesIndep = window[i].N[ig]/(window[i].tau[ig]/window[i].dt);
- window[i].sigma[ig] = sig * nSamplesIndep/(nSamplesIndep-1);
+ nSamplesIndep = window[i].N[ig] / (window[i].tau[ig] / window[i].dt);
+ window[i].sigma[ig] = sig * nSamplesIndep / (nSamplesIndep - 1);
}
else
{
/*! \brief
* Use histograms to compute average force on pull group.
*/
-static void computeAverageForce(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt)
+static void computeAverageForce(t_UmbrellaWindow* window, int nWindows, t_UmbrellaOptions* opt)
{
int i, j, bins = opt->bins, k;
double dz, min = opt->min, max = opt->max, displAv, temp, distance, ztot, ztot_half, w, weight;
double posmirrored;
- dz = (max-min)/bins;
- ztot = opt->max-min;
- ztot_half = ztot/2;
+ dz = (max - min) / bins;
+ ztot = opt->max - min;
+ ztot_half = ztot / 2;
/* Compute average displacement from histograms */
for (j = 0; j < nWindows; ++j)
weight = 0.0;
for (i = 0; i < opt->bins; ++i)
{
- temp = (1.0*i+0.5)*dz+min;
+ temp = (1.0 * i + 0.5) * dz + min;
distance = temp - window[j].pos[k];
if (opt->bCycl)
- { /* in cyclic wham: */
- if (distance > ztot_half) /* |distance| < ztot_half */
+ { /* in cyclic wham: */
+ if (distance > ztot_half) /* |distance| < ztot_half */
{
distance -= ztot;
}
distance += ztot;
}
}
- w = window[j].Histo[k][i]/window[j].g[k];
- displAv += w*distance;
- weight += w;
+ w = window[j].Histo[k][i] / window[j].g[k];
+ displAv += w * distance;
+ weight += w;
/* Are we near min or max? We are getting wrong forces from the histgrams since
the histograms are zero outside [min,max). Therefore, assume that the position
on the other side of the histomgram center is equally likely. */
if (!opt->bCycl)
{
- posmirrored = window[j].pos[k]-distance;
+ posmirrored = window[j].pos[k] - distance;
if (posmirrored >= max || posmirrored < min)
{
- displAv += -w*distance;
- weight += w;
+ displAv += -w * distance;
+ weight += w;
}
}
}
- displAv /= weight;
+ displAv /= weight;
/* average force from average displacement */
- window[j].forceAv[k] = displAv*window[j].k[k];
+ window[j].forceAv[k] = displAv * window[j].k[k];
/* sigma from average square displacement */
/* window[j].sigma [k] = sqrt(displAv2); */
/* printf("Win %d, sigma = %f\n",j,sqrt(displAv2)); */
/*! \brief
* Check if the complete reaction coordinate is covered by the histograms
*/
-static void checkReactionCoordinateCovered(t_UmbrellaWindow *window, int nwins,
- t_UmbrellaOptions *opt)
+static void checkReactionCoordinateCovered(t_UmbrellaWindow* window, int nwins, t_UmbrellaOptions* opt)
{
int i, ig, j, bins = opt->bins;
bool bBoundary;
count[j] += window[i].Histo[ig][j];
}
}
- avcount += 1.0*count[j];
+ avcount += 1.0 * count[j];
}
avcount /= bins;
for (j = 0; j < bins; ++j)
{
- relcount = count[j]/avcount;
- z = (j+0.5)*opt->dz+opt->min;
- bBoundary = j<bins/20 || (bins-j)>bins/20;
+ relcount = count[j] / avcount;
+ z = (j + 0.5) * opt->dz + opt->min;
+ bBoundary = j < bins / 20 || (bins - j) > bins / 20;
/* check for bins with no data */
if (count[j] == 0)
{
- fprintf(stderr, "\nWARNING, no data point in bin %d (z=%g) !\n"
- "You may not get a reasonable profile. Check your histograms!\n", j, z);
+ fprintf(stderr,
+ "\nWARNING, no data point in bin %d (z=%g) !\n"
+ "You may not get a reasonable profile. Check your histograms!\n",
+ j, z);
}
/* and check for poor sampling */
else if (relcount < 0.005 && !bBoundary)
*
* This speeds up the convergence by roughly a factor of 2
*/
-static void guessPotByIntegration(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt, const char *xlabel)
+static void guessPotByIntegration(t_UmbrellaWindow* window, int nWindows, t_UmbrellaOptions* opt, const char* xlabel)
{
int i, j, ig, bins = opt->bins, nHist, winmin, groupmin;
double dz, min = opt->min, *pot, pos, hispos, dist, diff, fAv, distmin, *f;
- FILE *fp;
+ FILE* fp;
- dz = (opt->max-min)/bins;
+ dz = (opt->max - min) / bins;
printf("Getting initial potential by integration.\n");
snew(f, bins);
for (j = 0; j < opt->bins; ++j)
{
- pos = (1.0*j+0.5)*dz+min;
+ pos = (1.0 * j + 0.5) * dz + min;
nHist = 0;
fAv = 0.;
distmin = 1e20;
for (ig = 0; ig < window[i].nPull; ig++)
{
hispos = window[i].pos[ig];
- dist = std::abs(hispos-pos);
+ dist = std::abs(hispos - pos);
/* average force within bin */
- if (dist < dz/2)
+ if (dist < dz / 2)
{
nHist++;
fAv += window[i].forceAv[ig];
/* if no histogram found in this bin, use closest histogram */
if (nHist > 0)
{
- fAv = fAv/nHist;
+ fAv = fAv / nHist;
}
else
{
}
for (j = 1; j < opt->bins; ++j)
{
- pot[j] = pot[j-1] - 0.5*dz*(f[j-1]+f[j]);
+ pot[j] = pot[j - 1] - 0.5 * dz * (f[j - 1] + f[j]);
}
/* cyclic wham: linearly correct possible offset */
if (opt->bCycl)
{
- diff = (pot[bins-1]-pot[0])/(bins-1);
+ diff = (pot[bins - 1] - pot[0]) / (bins - 1);
for (j = 1; j < opt->bins; ++j)
{
- pot[j] -= j*diff;
+ pot[j] -= j * diff;
}
}
if (opt->verbose)
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]);
+ fprintf(fp, "%g %g\n", (j + 0.5) * dz + opt->min, pot[j]);
}
xvgrclose(fp);
printf("verbose mode: wrote %s with PMF from interated forces\n", "pmfintegrated.xvg");
*/
for (j = 0; j < opt->bins; ++j)
{
- pot[j] = std::exp(-pot[j]/(BOLTZ*opt->Temperature));
+ pot[j] = std::exp(-pot[j] / (BOLTZ * opt->Temperature));
}
calc_z(pot, window, nWindows, opt, TRUE);
}
//! Count number of words in an line
-static int wordcount(char *ptr)
+static int wordcount(char* ptr)
{
int i, n, is[2];
int cur = 0;
for (i = 0; (ptr[i] != '\0'); i++)
{
is[cur] = isspace(ptr[i]);
- if ((i > 0) && (is[cur] && !is[1-cur]))
+ if ((i > 0) && (is[cur] && !is[1 - cur]))
{
n++;
}
- cur = 1-cur;
+ cur = 1 - cur;
}
return n;
}
*
* TO DO: ptr=fgets(...) is never freed (small memory leak)
*/
-static void readPullCoordSelection(t_UmbrellaOptions *opt, char **fnTpr, int nTpr)
+static void readPullCoordSelection(t_UmbrellaOptions* opt, char** fnTpr, int nTpr)
{
- FILE *fp;
+ FILE* fp;
int i, iline, n, len = STRLEN, temp;
char *ptr = nullptr, *tmpbuf = nullptr;
char fmt[1024], fmtign[1024];
snew(tmpbuf, len);
sizenow = 0;
iline = 0;
- while ( (ptr = fgets3(fp, tmpbuf, &len)) != nullptr)
+ while ((ptr = fgets3(fp, tmpbuf, &len)) != nullptr)
{
trim(ptr);
n = wordcount(ptr);
opt->nCoordsel = iline;
if (nTpr != opt->nCoordsel)
{
- gmx_fatal(FARGS, "Found %d tpr files but %d lines in %s\n", nTpr, opt->nCoordsel,
- opt->fnCoordSel);
+ gmx_fatal(FARGS, "Found %d tpr files but %d lines in %s\n", nTpr, opt->nCoordsel, opt->fnCoordSel);
}
printf("\nUse only these pull coordinates:\n");
for (iline = 0; iline < nTpr; iline++)
{
- printf("%s (%d of %d coordinates):", fnTpr[iline], opt->coordsel[iline].nUse, opt->coordsel[iline].n);
+ printf("%s (%d of %d coordinates):", fnTpr[iline], opt->coordsel[iline].nUse,
+ opt->coordsel[iline].n);
for (i = 0; i < opt->coordsel[iline].n; i++)
{
if (opt->coordsel[iline].bUse[i])
{
- printf(" %d", i+1);
+ printf(" %d", i + 1);
}
}
printf("\n");
}
//! Boolean XOR
-#define WHAMBOOLXOR(a, b) ( ((!(a)) && (b)) || ((a) && (!(b))))
+#define WHAMBOOLXOR(a, b) (((!(a)) && (b)) || ((a) && (!(b))))
//! Number of elements in fnm (used for command line parsing)
#define NFILE asize(fnm)
//! The main gmx wham routine
-int gmx_wham(int argc, char *argv[])
+int gmx_wham(int argc, char* argv[])
{
- const char *desc[] = {
+ const char* desc[] = {
"[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 ",
" provides the pull force output file names ([TT]pullf.xvg[tt]) with option [TT]-if[tt].",
" From the pull force the position in the umbrella potential is",
" computed. This does not work with tabulated umbrella potentials.",
- "* With option [TT]-ip[tt], the user provides file names of (gzipped) [REF].pdo[ref] files, i.e.",
+ "* With option [TT]-ip[tt], the user provides file names of (gzipped) [REF].pdo[ref] ",
+ " files, i.e.",
" the GROMACS 3.3 umbrella output files. If you have some unusual",
" reaction coordinate you may also generate your own [REF].pdo[ref] files and",
- " feed them with the [TT]-ip[tt] option into to [THISMODULE]. The [REF].pdo[ref] file header",
- " must be similar to the following::",
+ " feed them with the [TT]-ip[tt] option into to [THISMODULE]. The [REF].pdo[ref] file ",
+ " header must be similar to the following::",
"",
" # UMBRELLA 3.0",
" # Component selection: 0 0 1",
" a data column for each pull group follows (i.e. the displacement",
" with respect to the umbrella center). Up to four pull groups are possible ",
" per [REF].pdo[ref] file at present.[PAR]",
- "By default, all pull coordinates found in all pullx/pullf files are used in WHAM. If only ",
- "some of the pull coordinates should be used, a pull coordinate selection file (option [TT]-is[tt]) can ",
+ "By default, all pull coordinates found in all pullx/pullf files are used in WHAM. If ",
+ "only ",
+ "some of the pull coordinates should be used, a pull coordinate selection file (option ",
+ "[TT]-is[tt]) can ",
"be provided. The selection file must contain one line for each tpr file in tpr-files.dat.",
- "Each of these lines must contain one digit (0 or 1) for each pull coordinate in the tpr file. ",
- "Here, 1 indicates that the pull coordinate is used in WHAM, and 0 means it is omitted. Example:",
- "If you have three tpr files, each containing 4 pull coordinates, but only pull coordinates 1 and 2 should be ",
+ "Each of these lines must contain one digit (0 or 1) for each pull coordinate in the tpr ",
+ "file. ",
+ "Here, 1 indicates that the pull coordinate is used in WHAM, and 0 means it is omitted. ",
+ "Example:",
+ "If you have three tpr files, each containing 4 pull coordinates, but only pull ",
+ "coordinates 1 and 2 should be ",
"used, coordsel.dat looks like this::",
"",
" 1 1 0 0",
"",
"Always check whether the histograms sufficiently overlap.[PAR]",
"The umbrella potential is assumed to be harmonic and the force constants are ",
- "read from the [REF].tpr[ref] or [REF].pdo[ref] files. If a non-harmonic umbrella force was applied ",
+ "read from the [REF].tpr[ref] or [REF].pdo[ref] files. If a non-harmonic umbrella force ",
+ "was applied ",
"a tabulated potential can be provided with [TT]-tab[tt].",
"",
"WHAM options",
"less robust) method such as fitting to a double exponential, you can ",
"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 ([REF].pdo[ref] or pullx/f file) and one column per pull coordinate in the respective file.",
+ "input file ([REF].pdo[ref] or pullx/f file) and one column per pull coordinate in the ",
+ "respective file.",
"",
"Error analysis",
"^^^^^^^^^^^^^^",
" and properly autocorrelated. The autocorrelation time (ACT) for each window must be ",
" known, so use [TT]-ac[tt] or provide the ACT with [TT]-iiact[tt]. If the ACT of all ",
" windows are identical (and known), you can also provide them with [TT]-bs-tau[tt]. ",
- " Note that this method may severely underestimate the error in case of limited sampling, ",
+ " Note that this method may severely underestimate the error in case of limited ",
+ " sampling, ",
" that is if individual histograms do not represent the complete phase space at ",
" the respective positions.",
"* [TT]traj-gauss[tt] The same as method [TT]traj[tt], but the trajectories are ",
"the histograms."
};
- const char *en_unit[] = {nullptr, "kJ", "kCal", "kT", nullptr};
- const char *en_unit_label[] = {"", "E (kJ mol\\S-1\\N)", "E (kcal mol\\S-1\\N)", "E (kT)", nullptr};
- const char *en_bsMethod[] = { nullptr, "b-hist", "hist", "traj", "traj-gauss", nullptr };
+ const char* en_unit[] = { nullptr, "kJ", "kCal", "kT", nullptr };
+ const char* en_unit_label[] = { "", "E (kJ mol\\S-1\\N)", "E (kcal mol\\S-1\\N)", "E (kT)", nullptr };
+ const char* en_bsMethod[] = { nullptr, "b-hist", "hist", "traj", "traj-gauss", nullptr };
static t_UmbrellaOptions opt;
- t_pargs pa[] = {
- { "-min", FALSE, etREAL, {&opt.min},
- "Minimum coordinate in profile"},
- { "-max", FALSE, etREAL, {&opt.max},
- "Maximum coordinate in profile"},
- { "-auto", FALSE, etBOOL, {&opt.bAuto},
- "Determine min and max automatically"},
- { "-bins", FALSE, etINT, {&opt.bins},
- "Number of bins in profile"},
- { "-temp", FALSE, etREAL, {&opt.Temperature},
- "Temperature"},
- { "-tol", FALSE, etREAL, {&opt.Tolerance},
- "Tolerance"},
- { "-v", FALSE, etBOOL, {&opt.verbose},
- "Verbose mode"},
- { "-b", FALSE, etREAL, {&opt.tmin},
- "First time to analyse (ps)"},
- { "-e", FALSE, etREAL, {&opt.tmax},
- "Last time to analyse (ps)"},
- { "-dt", FALSE, etREAL, {&opt.dt},
- "Analyse only every dt ps"},
- { "-histonly", FALSE, etBOOL, {&opt.bHistOnly},
- "Write histograms and exit"},
- { "-boundsonly", FALSE, etBOOL, {&opt.bBoundsOnly},
- "Determine min and max and exit (with [TT]-auto[tt])"},
- { "-log", FALSE, etBOOL, {&opt.bLog},
- "Calculate the log of the profile before printing"},
- { "-unit", FALSE, etENUM, {en_unit},
- "Energy unit in case of log output" },
- { "-zprof0", FALSE, etREAL, {&opt.zProf0},
- "Define profile to 0.0 at this position (with [TT]-log[tt])"},
- { "-cycl", FALSE, etBOOL, {&opt.bCycl},
- "Create cyclic/periodic profile. Assumes min and max are the same point."},
- { "-sym", FALSE, etBOOL, {&opt.bSym},
- "Symmetrize profile around z=0"},
- { "-hist-eq", FALSE, etBOOL, {&opt.bHistEq},
- "HIDDENEnforce equal weight for all histograms. (Non-Weighed-HAM)"},
- { "-ac", FALSE, etBOOL, {&opt.bCalcTauInt},
- "Calculate integrated autocorrelation times and use in wham"},
- { "-acsig", FALSE, etREAL, {&opt.sigSmoothIact},
- "Smooth autocorrelation times along reaction coordinate with Gaussian of this [GRK]sigma[grk]"},
- { "-ac-trestart", FALSE, etREAL, {&opt.acTrestart},
- "When computing autocorrelation functions, restart computing every .. (ps)"},
- { "-acred", FALSE, etBOOL, {&opt.bAllowReduceIact},
+ t_pargs pa[] = {
+ { "-min", FALSE, etREAL, { &opt.min }, "Minimum coordinate in profile" },
+ { "-max", FALSE, etREAL, { &opt.max }, "Maximum coordinate in profile" },
+ { "-auto", FALSE, etBOOL, { &opt.bAuto }, "Determine min and max automatically" },
+ { "-bins", FALSE, etINT, { &opt.bins }, "Number of bins in profile" },
+ { "-temp", FALSE, etREAL, { &opt.Temperature }, "Temperature" },
+ { "-tol", FALSE, etREAL, { &opt.Tolerance }, "Tolerance" },
+ { "-v", FALSE, etBOOL, { &opt.verbose }, "Verbose mode" },
+ { "-b", FALSE, etREAL, { &opt.tmin }, "First time to analyse (ps)" },
+ { "-e", FALSE, etREAL, { &opt.tmax }, "Last time to analyse (ps)" },
+ { "-dt", FALSE, etREAL, { &opt.dt }, "Analyse only every dt ps" },
+ { "-histonly", FALSE, etBOOL, { &opt.bHistOnly }, "Write histograms and exit" },
+ { "-boundsonly",
+ FALSE,
+ etBOOL,
+ { &opt.bBoundsOnly },
+ "Determine min and max and exit (with [TT]-auto[tt])" },
+ { "-log", FALSE, etBOOL, { &opt.bLog }, "Calculate the log of the profile before printing" },
+ { "-unit", FALSE, etENUM, { en_unit }, "Energy unit in case of log output" },
+ { "-zprof0",
+ FALSE,
+ etREAL,
+ { &opt.zProf0 },
+ "Define profile to 0.0 at this position (with [TT]-log[tt])" },
+ { "-cycl",
+ FALSE,
+ etBOOL,
+ { &opt.bCycl },
+ "Create cyclic/periodic profile. Assumes min and max are the same point." },
+ { "-sym", FALSE, etBOOL, { &opt.bSym }, "Symmetrize profile around z=0" },
+ { "-hist-eq",
+ FALSE,
+ etBOOL,
+ { &opt.bHistEq },
+ "HIDDENEnforce equal weight for all histograms. (Non-Weighed-HAM)" },
+ { "-ac",
+ FALSE,
+ etBOOL,
+ { &opt.bCalcTauInt },
+ "Calculate integrated autocorrelation times and use in wham" },
+ { "-acsig",
+ FALSE,
+ etREAL,
+ { &opt.sigSmoothIact },
+ "Smooth autocorrelation times along reaction coordinate with Gaussian of this "
+ "[GRK]sigma[grk]" },
+ { "-ac-trestart",
+ FALSE,
+ etREAL,
+ { &opt.acTrestart },
+ "When computing autocorrelation functions, restart computing every .. (ps)" },
+ { "-acred",
+ FALSE,
+ etBOOL,
+ { &opt.bAllowReduceIact },
"HIDDENWhen smoothing the ACTs, allows one to reduce ACTs. Otherwise, only increase ACTs "
- "during smoothing"},
- { "-nBootstrap", FALSE, etINT, {&opt.nBootStrap},
+ "during smoothing" },
+ { "-nBootstrap",
+ FALSE,
+ etINT,
+ { &opt.nBootStrap },
"nr of bootstraps to estimate statistical uncertainty (e.g., 200)" },
- { "-bs-method", FALSE, etENUM, {en_bsMethod},
- "Bootstrap method" },
- { "-bs-tau", FALSE, etREAL, {&opt.tauBootStrap},
- "Autocorrelation time (ACT) assumed for all histograms. Use option [TT]-ac[tt] if ACT is unknown."},
- { "-bs-seed", FALSE, etINT, {&opt.bsSeed},
- "Seed for bootstrapping. (-1 = use time)"},
- { "-histbs-block", FALSE, etINT, {&opt.histBootStrapBlockLength},
- "When mixing histograms only mix within blocks of [TT]-histbs-block[tt]."},
- { "-vbs", FALSE, etBOOL, {&opt.bs_verbose},
- "Verbose bootstrapping. Print the CDFs and a histogram file for each bootstrap."},
- { "-stepout", FALSE, etINT, {&opt.stepchange},
- "HIDDENWrite maximum change every ... (set to 1 with [TT]-v[tt])"},
- { "-updateContr", FALSE, etINT, {&opt.stepUpdateContrib},
- "HIDDENUpdate table with significan contributions to WHAM every ... iterations"},
+ { "-bs-method", FALSE, etENUM, { en_bsMethod }, "Bootstrap method" },
+ { "-bs-tau",
+ FALSE,
+ etREAL,
+ { &opt.tauBootStrap },
+ "Autocorrelation time (ACT) assumed for all histograms. Use option [TT]-ac[tt] if ACT is "
+ "unknown." },
+ { "-bs-seed", FALSE, etINT, { &opt.bsSeed }, "Seed for bootstrapping. (-1 = use time)" },
+ { "-histbs-block",
+ FALSE,
+ etINT,
+ { &opt.histBootStrapBlockLength },
+ "When mixing histograms only mix within blocks of [TT]-histbs-block[tt]." },
+ { "-vbs",
+ FALSE,
+ etBOOL,
+ { &opt.bs_verbose },
+ "Verbose bootstrapping. Print the CDFs and a histogram file for each bootstrap." },
+ { "-stepout",
+ FALSE,
+ etINT,
+ { &opt.stepchange },
+ "HIDDENWrite maximum change every ... (set to 1 with [TT]-v[tt])" },
+ { "-updateContr",
+ FALSE,
+ etINT,
+ { &opt.stepUpdateContrib },
+ "HIDDENUpdate table with significan contributions to WHAM every ... iterations" },
};
- t_filenm fnm[] = {
- { efDAT, "-ix", "pullx-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs */
- { efDAT, "-if", "pullf-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs */
- { efDAT, "-it", "tpr-files", ffOPTRD}, /* wham input: tprs */
- { efDAT, "-ip", "pdo-files", ffOPTRD}, /* wham input: pdo files (gmx3 style) */
- { efDAT, "-is", "coordsel", ffOPTRD}, /* input: select pull coords to use */
- { efXVG, "-o", "profile", ffWRITE }, /* output file for profile */
- { efXVG, "-hist", "histo", ffWRITE}, /* output file for histograms */
- { efXVG, "-oiact", "iact", ffOPTWR}, /* writing integrated autocorrelation times */
- { efDAT, "-iiact", "iact-in", ffOPTRD}, /* reading integrated autocorrelation times */
- { efXVG, "-bsres", "bsResult", ffOPTWR}, /* average and errors of bootstrap analysis */
- { efXVG, "-bsprof", "bsProfs", ffOPTWR}, /* output file for bootstrap profiles */
- { efDAT, "-tab", "umb-pot", ffOPTRD}, /* Tabulated umbrella potential (if not harmonic) */
+ t_filenm fnm[] = {
+ { efDAT, "-ix", "pullx-files", ffOPTRD }, /* wham input: pullf.xvg's and tprs */
+ { efDAT, "-if", "pullf-files", ffOPTRD }, /* wham input: pullf.xvg's and tprs */
+ { efDAT, "-it", "tpr-files", ffOPTRD }, /* wham input: tprs */
+ { efDAT, "-ip", "pdo-files", ffOPTRD }, /* wham input: pdo files (gmx3 style) */
+ { efDAT, "-is", "coordsel", ffOPTRD }, /* input: select pull coords to use */
+ { efXVG, "-o", "profile", ffWRITE }, /* output file for profile */
+ { efXVG, "-hist", "histo", ffWRITE }, /* output file for histograms */
+ { efXVG, "-oiact", "iact", ffOPTWR }, /* writing integrated autocorrelation times */
+ { efDAT, "-iiact", "iact-in", ffOPTRD }, /* reading integrated autocorrelation times */
+ { efXVG, "-bsres", "bsResult", ffOPTWR }, /* average and errors of bootstrap analysis */
+ { efXVG, "-bsprof", "bsProfs", ffOPTWR }, /* output file for bootstrap profiles */
+ { efDAT, "-tab", "umb-pot", ffOPTRD }, /* Tabulated umbrella potential (if not harmonic) */
};
- int i, j, l, nfiles, nwins, nfiles2;
- t_UmbrellaHeader header;
- t_UmbrellaWindow * window = nullptr;
- double *profile, maxchange = 1e20;
- gmx_bool bMinSet, bMaxSet, bAutoSet, bExact = FALSE;
- char **fninTpr, **fninPull, **fninPdo;
- const char *fnPull;
- FILE *histout, *profout;
- char xlabel[STRLEN], ylabel[256], title[256];
+ int i, j, l, nfiles, nwins, nfiles2;
+ t_UmbrellaHeader header;
+ t_UmbrellaWindow* window = nullptr;
+ double * profile, maxchange = 1e20;
+ gmx_bool bMinSet, bMaxSet, bAutoSet, bExact = FALSE;
+ char ** fninTpr, **fninPull, **fninPdo;
+ const char* fnPull;
+ FILE * histout, *profout;
+ char xlabel[STRLEN], ylabel[256], title[256];
opt.bins = 200;
opt.verbose = FALSE;
opt.stepchange = 100;
opt.stepUpdateContrib = 100;
- if (!parse_common_args(&argc, argv, 0,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &opt.oenv))
+ if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr,
+ &opt.oenv))
{
return 0;
}
opt.unit = nenum(en_unit);
opt.bsMethod = nenum(en_bsMethod);
- opt.bProf0Set = opt2parg_bSet("-zprof0", asize(pa), pa);
+ opt.bProf0Set = opt2parg_bSet("-zprof0", asize(pa), pa);
opt.bTab = opt2bSet("-tab", NFILE, fnm);
opt.bPdo = opt2bSet("-ip", NFILE, fnm);
opt.bPullx = opt2bSet("-ix", NFILE, fnm);
opt.bPullf = opt2bSet("-if", NFILE, fnm);
opt.bTauIntGiven = opt2bSet("-iiact", NFILE, fnm);
- if (opt.bTab && opt.bPullf)
+ if (opt.bTab && opt.bPullf)
{
- gmx_fatal(FARGS, "Force input does not work with tabulated potentials. "
+ gmx_fatal(FARGS,
+ "Force input does not work with tabulated potentials. "
"Provide pullx.xvg or pdo files!");
}
}
if (!opt.bPdo && !(opt.bTpr || opt.bPullf || opt.bPullx))
{
- gmx_fatal(FARGS, "gmx wham supports three input modes, pullx, pullf, or pdo file input."
+ gmx_fatal(FARGS,
+ "gmx wham supports three input modes, pullx, pullf, or pdo file input."
"\n\n Check gmx wham -h !");
}
opt.fnPullx = opt2fn("-ix", NFILE, fnm);
opt.fnCoordSel = opt2fn_null("-is", NFILE, fnm);
- bMinSet = opt2parg_bSet("-min", asize(pa), pa);
- bMaxSet = opt2parg_bSet("-max", asize(pa), pa);
- bAutoSet = opt2parg_bSet("-auto", asize(pa), pa);
- if ( (bMinSet || bMaxSet) && bAutoSet)
+ bMinSet = opt2parg_bSet("-min", asize(pa), pa);
+ bMaxSet = opt2parg_bSet("-max", asize(pa), pa);
+ bAutoSet = opt2parg_bSet("-auto", asize(pa), pa);
+ if ((bMinSet || bMaxSet) && bAutoSet)
{
gmx_fatal(FARGS, "With -auto, do not give -min or -max\n");
}
- if ( (bMinSet && !bMaxSet) || (!bMinSet && bMaxSet))
+ if ((bMinSet && !bMaxSet) || (!bMinSet && bMaxSet))
{
gmx_fatal(FARGS, "When giving -min, you must give -max (and vice versa), too\n");
}
if (opt.bTauIntGiven && opt.bCalcTauInt)
{
- gmx_fatal(FARGS, "Either read (option -iiact) or calculate (option -ac) the\n"
+ gmx_fatal(FARGS,
+ "Either read (option -iiact) or calculate (option -ac) the\n"
"the autocorrelation times. Not both.");
}
if (opt.tauBootStrap > 0.0 && opt2parg_bSet("-ac", asize(pa), pa))
{
- gmx_fatal(FARGS, "Either compute autocorrelation times (ACTs) (option -ac) or "
+ gmx_fatal(FARGS,
+ "Either compute autocorrelation times (ACTs) (option -ac) or "
"provide it with -bs-tau for bootstrapping. Not Both.\n");
}
if (opt.tauBootStrap > 0.0 && opt2bSet("-iiact", NFILE, fnm))
{
- gmx_fatal(FARGS, "Either provide autocorrelation times (ACTs) with file iact-in.dat "
+ gmx_fatal(FARGS,
+ "Either provide autocorrelation times (ACTs) with file iact-in.dat "
"(option -iiact) or define all ACTs with -bs-tau for bootstrapping\n. Not Both.");
}
fnPull = opt.bPullf ? opt.fnPullf : opt.fnPullx;
read_wham_in(fnPull, &fninPull, &nfiles2, &opt);
- printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n",
- nfiles, nfiles2, opt.bPullf ? "force" : "position", opt.fnTpr, fnPull);
+ printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n", nfiles, nfiles2,
+ opt.bPullf ? "force" : "position", opt.fnTpr, fnPull);
if (nfiles != nfiles2)
{
- gmx_fatal(FARGS, "Found %d file names in %s, but %d in %s\n", nfiles,
- opt.fnTpr, nfiles2, fnPull);
+ gmx_fatal(FARGS, "Found %d file names in %s, but %d in %s\n", nfiles, opt.fnTpr, nfiles2, fnPull);
}
/* Read file that selects the pull group to be used */
read_tpr_pullxf_files(fninTpr, fninPull, nfiles, &header, window, &opt);
}
else
- { /* reading pdo files */
- if (opt.fnCoordSel != nullptr)
+ { /* reading pdo files */
+ if (opt.fnCoordSel != nullptr)
{
- gmx_fatal(FARGS, "Reading a -is file is not supported with PDO input files.\n"
- "Use awk or a similar tool to pick the required pull groups from your PDO files\n");
+ gmx_fatal(FARGS,
+ "Reading a -is file is not supported with PDO input files.\n"
+ "Use awk or a similar tool to pick the required pull groups from your PDO "
+ "files\n");
}
read_wham_in(opt.fnPdo, &fninPdo, &nfiles, &opt);
printf("Found %d pdo files in %s\n", nfiles, opt.fnPdo);
}
/* write histograms */
- histout = xvgropen(opt2fn("-hist", NFILE, fnm), "Umbrella histograms",
- xlabel, "count", opt.oenv);
+ histout = xvgropen(opt2fn("-hist", NFILE, fnm), "Umbrella histograms", xlabel, "count", opt.oenv);
for (l = 0; l < opt.bins; ++l)
{
- fprintf(histout, "%e\t", (l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
+ fprintf(histout, "%e\t", (l + 0.5) / opt.bins * (opt.max - opt.min) + opt.min);
for (i = 0; i < nwins; ++i)
{
for (j = 0; j < window[i].nPull; ++j)
i = 0;
do
{
- if ( (i%opt.stepUpdateContrib) == 0)
+ if ((i % opt.stepUpdateContrib) == 0)
{
setup_acc_wham(profile, window, nwins, &opt);
}
printf("Switched to exact iteration in iteration %d\n", i);
}
calc_profile(profile, window, nwins, &opt, bExact);
- if (((i%opt.stepchange) == 0 || i == 1) && i != 0)
+ if (((i % opt.stepchange) == 0 || i == 1) && i != 0)
{
printf("\t%4d) Maximum change %e\n", i, maxchange);
}
i++;
- }
- while ( (maxchange = calc_z(profile, window, nwins, &opt, bExact)) > opt.Tolerance || !bExact);
+ } while ((maxchange = calc_z(profile, window, nwins, &opt, bExact)) > opt.Tolerance || !bExact);
printf("Converged in %d iterations. Final maximum change %g\n", i, maxchange);
/* calc error from Kumar's formula */
profout = xvgropen(opt2fn("-o", NFILE, fnm), title, xlabel, ylabel, opt.oenv);
for (i = 0; i < opt.bins; ++i)
{
- fprintf(profout, "%e\t%e\n", (i+0.5)/opt.bins*(opt.max-opt.min)+opt.min, profile[i]);
+ fprintf(profout, "%e\t%e\n", (i + 0.5) / opt.bins * (opt.max - opt.min) + opt.min, profile[i]);
}
xvgrclose(profout);
printf("Wrote %s\n", opt2fn("-o", NFILE, fnm));
if (opt.nBootStrap)
{
do_bootstrapping(opt2fn("-bsres", NFILE, fnm), opt2fn("-bsprof", NFILE, fnm),
- opt2fn("-hist", NFILE, fnm),
- xlabel, ylabel, profile, window, nwins, &opt);
+ opt2fn("-hist", NFILE, fnm), xlabel, ylabel, profile, window, nwins, &opt);
}
sfree(profile);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff