diff --git a/README.md b/README.md index d7fe6e0..b4c69a5 100644 --- a/README.md +++ b/README.md @@ -62,7 +62,7 @@ The first column will be ignored and is followed by N reaction coordinates x. Shipped under the GPLv3 license. USAGE: - wham [FLAGS] [OPTIONS] --bins --max --file --min --temperature + wham [FLAGS] [OPTIONS] --bins --max --file --min --temperature FLAGS: -c, --cyclic For periodic reaction coordinates. If this is set, the first and last coordinate bin in each @@ -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*). diff --git a/src/cli.yml b/src/cli.yml index 216e506..9b1f067 100644 --- a/src/cli.yml +++ b/src/cli.yml @@ -1,6 +1,6 @@ name: wham version: "1.0.0" -author: D. Bauer +author: D. Bauer 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. diff --git a/src/correlation_analysis.rs b/src/correlation_analysis.rs index 6ff24d0..fea772c 100644 --- a/src/correlation_analysis.rs +++ b/src/correlation_analysis.rs @@ -70,13 +70,13 @@ mod tests { let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg"); let autocorr = super::autocorrelation(×eries); 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); } diff --git a/src/io.rs b/src/io.rs index 309c015..2b8188e 100644 --- a/src/io.rs +++ b/src/io.rs @@ -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>> { Ok(timeseries) } + // calculates the inefficiency for every collective variable // filters the timeseries based on the highest inefficiency fn uncorrelate(timeseries: Vec>, cfg: &Config) -> Vec> {