bsMethod_trajGauss
};
-//! Parameters of one pull coodinate
+//! Parameters of one pull coordinate
typedef struct
{
PullingAlgorithm pull_type; //!< such as constraint, umbrella, ...
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* 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.
*
* \name Basic WHAM options
*/
/*!\{*/
- 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
+ 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
}
}
/* Note: there are two contributions to bin k in the wham equations:
- i) N[j]*exp(- U/(BOLTZ*opt->Temperature) + window[i].z[j])
- ii) exp(- U/(BOLTZ*opt->Temperature))
+ i) N[j]*exp(- U/(c_boltz*opt->Temperature) + window[i].z[j])
+ ii) exp(- U/(c_boltz*opt->Temperature))
where U is the umbrella potential
If any of these number is larger wham_contrib_lim, I set contrib=TRUE
*/
{
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 / (gmx::c_boltz * opt->Temperature));
+ contrib2 = window[i].N[j]
+ * std::exp(-U / (gmx::c_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])
U = tabulated_pot(distance, opt); /* Use tabulated potential */
}
denom += invg * window[j].N[k]
- * std::exp(-U / (BOLTZ * opt->Temperature) + window[j].z[k]);
+ * std::exp(-U / (gmx::c_boltz * opt->Temperature) + window[j].z[k]);
}
}
profile[i] = num / denom;
{
U = tabulated_pot(distance, opt); /* Use tabulated potential */
}
- total += profile[k] * std::exp(-U / (BOLTZ * opt->Temperature));
+ total += profile[k] * std::exp(-U / (gmx::c_boltz * opt->Temperature));
}
/* Avoid floating point exception if window is far outside min and max */
if (total != 0.0)
}
else if (opt->unit == en_kJ)
{
- unit_factor = BOLTZ * opt->Temperature;
+ unit_factor = gmx::c_boltz * opt->Temperature;
}
else if (opt->unit == en_kCal)
{
- unit_factor = (BOLTZ / CAL2JOULE) * opt->Temperature;
+ unit_factor = (gmx::c_boltz / gmx::c_cal2Joule) * opt->Temperature;
}
else
{
* 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].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];
- 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 != FreeEnergyPerturbationType::No && ir->pull->coord[i].k != ir->pull->coord[i].kB)
{
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 "
+ "If you have umbrella and non-umbrella coordinates, you can select the "
"umbrella coordinates with gmx wham -is\n",
fn,
i + 1,
}
if (bFirst || opt->verbose)
{
+ // Note the comment above about the pull file column format!
printf("\nReading pull %s file %s, expecting %d columns:\n", quantity, fn, nColExpect);
+ printf("\tColumn for time: 1\n");
+ int nextColumn = 2;
for (i = 0; i < header->npullcrds; i++)
{
- 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",
- 1,
- nColCOMCrd[i],
- (header->bPrintRefValue ? "Yes" : "No"));
+ printf("\tColumn(s) with data for pull coordinate %d are\n", i + 1);
+ if (nColCOMCrd[i] > 0)
+ {
+ const int firstColumnForCOM = nextColumn;
+ const int lastColumnForCOM = firstColumnForCOM + nColCOMCrd[i] - 1;
+ const int columnForThisCoordinate = lastColumnForCOM + 1;
+ printf("\t\treaction coordinate: %d\n"
+ "\t\tcenter-of-mass of groups: %d through %d\n",
+ columnForThisCoordinate,
+ firstColumnForCOM,
+ lastColumnForCOM);
+ nextColumn = columnForThisCoordinate + 1;
+ }
+ else
+ {
+ const int columnForThisCoordinate = nextColumn;
+ printf("\t\treaction coordinate: %d\n"
+ "\t\tcenter-of-mass of groups: No\n",
+ columnForThisCoordinate);
+ ++nextColumn;
+ }
+ if (header->bPrintRefValue)
+ {
+ const int columnForRefValue = nextColumn;
+ printf("\t\treference position column: %d\n", columnForRefValue);
+ ++nextColumn;
+ }
+ else
+ {
+ printf("\t\treference position column: No\n");
+ }
}
printf("\tFound %d times in %s\n", nt, fn);
bFirst = FALSE;
int i;
real mintmp, maxtmp;
- printf("Reading %d tpr and pullf files\n", nfiles);
+ printf("Reading %d tpr and pullx/pullf files\n", nfiles);
/* min and max not given? */
if (opt->bAuto)
*/
for (j = 0; j < opt->bins; ++j)
{
- pot[j] = std::exp(-pot[j] / (BOLTZ * opt->Temperature));
+ pot[j] = std::exp(-pot[j] / (gmx::c_boltz * opt->Temperature));
}
calc_z(pot, window, nWindows, opt, TRUE);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff