version up

This commit is contained in:
Daniel Bauer
2020-10-26 10:52:22 +01:00
4 changed files with 23 additions and 21 deletions

View File

@@ -90,7 +90,16 @@ OPTIONS:
tolerance (defaults to 0.000001).
```
To run the two dimensional example (simulation of dialanine phi and psi angle):
Examples
---
The example folder contains input and output files for two simple test systems:
- 1d_cyclic: Phi torsion angle of dialanine in vaccum
- 2d_cyclic: Phi and psi torsion angles of the same system
The command below will run the two dimensional example (simulation of dialanine phi and psi angle) and calculate the free energy based on the two collective variables
in the range of -3.14 to 3.14, with 100 bins in each dimension and periodic collective variables:
```bash
wham --max 3.14,3.14 --min -3.14,-3.14 -T 300 --bins 100,100 --cyclic -f example/2d/metadata.dat
> Supplied WHAM options: Metadata=example/2d/metadata.dat, hist_min=[-3.14, -3.14], hist_max=[3.14, 3.14], bins=[100, 100] verbose=false, tolerance=0.000001, iterations=100000, temperature=300, cyclic=true
@@ -107,7 +116,6 @@ wham --max 3.14,3.14 --min -3.14,-3.14 -T 300 --bins 100,100 --cyclic -f example
```
After convergence, final bias offsets (F) and the free energy will be dumped to stdout and the output file is written.
The output file contains the free energy and probability for each bin. Probabilities are normalized to sum to P=1.0 and
the smallest free energy is set to 0 (with other free energies based on that).
```
@@ -127,7 +135,7 @@ the smallest free energy is set to 0 (with other free energies based on that).
Error analysis
---
WHAM can perform error analysis using the bayesian bootstrapping method. Every simulation window is assumed to be an
individual set of data point. By calculating probabilities N times with randomly assigned weights for each window,
individual set of data points. By calculating probabilities N times with randomly assigned weights for each window,
one can estimate the error as standard deviation between the N bootstrapping runs. For more details see
*Van der Spoel, D. et al. (2010). g_wham—A Free Weighted Histogram Analysis Implementation Including Robust Error and
Autocorrelation Estimates, JCTC, 6(12), 3713-3720*.
@@ -148,14 +156,6 @@ timeseries is used for unbiasing. A more detailed description of the method can
tempering simulations, JCTC 3(1):26-41*
Examples
---
The example folder contains input and output files for two simple test systems:
- 1d_cyclic: Phi torsion angle of dialanine in vaccum
- 2d_cyclic: Phi and psi torsion angles of the same system
TODO
---
- Replica exchange
@@ -165,7 +165,7 @@ License & Citing
WHAM is licensed under the GPL-3.0 license. Please read the LICENSE file in this
repository for more information.
There's no publication for this WHAM implementation. However, there is a citeabe DOI. If you use this software for your work, please consider citing it: *Bauer, D, WHAM - An efficient weighted histogram analysis implementation written in Rust, Zenodo. https://doi.org/10.5281/zenodo.1488597*
There's no publication for this WHAM implementation. However, there is a citeabe DOI. If you use this software for your work, please consider citing it: *Bauer, D., WHAM - An efficient weighted histogram analysis implementation written in Rust, Zenodo. https://doi.org/10.5281/zenodo.1488597*
Parts of this work, especially some perfomance optimizations and the I/O format, are inspired by the
implementation of A. Grossfield (*Grossfield, A, WHAM: the weighted histogram analysis method, http://membrane.urmc.rochester.edu/content/wham*).

View File

@@ -1,6 +1,6 @@
name: wham
version: "1.0.0"
author: D. Bauer <bauer@bio.tu-darmstadt.de>
author: D. Bauer <bauer@cbs.tu-darmstadt.de>
about: |
wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports potential of mean force (PMF) calculations in multiple dimensions at constant temperature.

View File

@@ -70,13 +70,13 @@ mod tests {
let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
let autocorr = super::autocorrelation(&timeseries);
let expected = [
0.6919008655979143, 0.5331719399671355,
0.20472620956463589, -0.002850876458920514,
-0.13842850077938146, -0.2652923552973232,
-0.31427198272235385, -0.2617505151557693,
-0.20594864730290338, -0.13310019091811812,
-0.1887568901193426, -0.1944936625424933,
-0.1963599673189461, -0.11391587833838003];
0.691_900_865_597_914_3, 0.533_171_939_967_135_5,
0.204_726_209_564_635_89, -0.002_850_876_458_920_514,
-0.138_428_500_779_381_46, -0.265_292_355_297_323_2,
-0.314_271_982_722_353_85, -0.261_750_515_155_769_3,
-0.205_948_647_302_903_38, -0.133_100_190_918_118_12,
-0.188_756_890_119_342_6, -0.194_493_662_542_493_3,
-0.196_359_967_318_946_1, -0.113_915_878_338_380_03];
for (actual, expected) in autocorr.iter().zip(expected.iter()) {
assert!((actual-expected).abs() < 0.001);
}

View File

@@ -132,6 +132,7 @@ fn is_in_time_boundaries(time: f64, cfg: &Config) -> bool {
false
}
// parse a time series file into a histogram
fn read_window_file(window_file: &str, cfg: &Config) -> Result<(Histogram, usize)> {
// total number of bins is the product of all dimensions length
@@ -210,6 +211,7 @@ fn read_timeseries(window_file: &str, cfg: &Config) -> Result<Vec<Vec<f64>>> {
Ok(timeseries)
}
// calculates the inefficiency for every collective variable
// filters the timeseries based on the highest inefficiency
fn uncorrelate(timeseries: Vec<Vec<f64>>, cfg: &Config) -> Vec<Vec<f64>> {