From 2cc34919b01b5cedbed570c2d9f493ef89598834 Mon Sep 17 00:00:00 2001 From: Daniel Bauer Date: Sun, 25 Oct 2020 22:40:13 +0100 Subject: [PATCH 1/3] Autocorrelation (#2) * calculate autocorrelation stats_ineff and tau * uncorr flag * read timeseries into vector * uncorrelate data * README * Update README.md * Update README.md Co-authored-by: Daniel Bauer --- README.md | 23 +++- example/1d_cyclic/wham_uncorrelated.out | 101 ++++++++++++++++ src/cli.yml | 6 + src/correlation_analysis.rs | 102 +++++++++++++++++ src/error_analysis.rs | 8 +- src/io.rs | 146 ++++++++++++++++++------ src/lib.rs | 6 +- src/main.rs | 4 +- tests/examples.rs | 20 ++++ 9 files changed, 370 insertions(+), 46 deletions(-) create mode 100644 example/1d_cyclic/wham_uncorrelated.out create mode 100644 src/correlation_analysis.rs diff --git a/README.md b/README.md index 6553ccd..3793af9 100644 --- a/README.md +++ b/README.md @@ -12,9 +12,10 @@ from umbrella sampling simulations. For more details on the method, I suggest *R Features --- - Fast, especially for small systems -- Multithreaded -- Multidimensional -- Error analysis +- Multithreaded (automatically runs on all available cores) +- Multidimensional (any number of collective variables are possible) +- Autocorrelation to remove correlated samples +- Error analysis via bootstrapping - Unit tested Installation @@ -67,6 +68,8 @@ FLAGS: -c, --cyclic For periodic reaction coordinates. If this is set, the first and last coordinate bin in each dimension are treated as neighbors for the bias calculation. -h, --help Prints help information + -g, --uncorr Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated + samples (default is off). -V, --version Prints version information -v, --verbose Enables verbose output. @@ -85,7 +88,6 @@ OPTIONS: -T, --temperature WHAM temperature in Kelvin. -t, --tolerance Abortion criteria for WHAM calculation. WHAM stops if abs(F_new - F_old) < tolerance (defaults to 0.000001). - ``` To run the two dimensional example (simulation of dialanine phi and psi angle): @@ -134,6 +136,18 @@ To perform bayesian bootstrapping in WHAM, use the ```-bt ``` flag to perf runs. The error estimates of bin probabilities and free energy will be given as standard error (SE) in a separate column (+/-) in the output file. If no error analysis is performed, these columns are set to 0.0. +Autocorrelation analysis +--- +With the ```--uncorr``` flag, WHAM calculates the autocorrelation time ```tau``` for all timeseries and all collective +variables. Timeseries are then filtered based on their highest autocorrelation time to remove correlated samples from +the dataset. This reduces the number of data points but can improve the accuracy of the result. + +For filtering, the statistical inefficiency `g` is calculated: ```g = 1 + 2*tau```, and only every `g`th element of the +timeseries is used for unbiasing. A more detailed description of the method can be found in +*Chodera, J.D. et al. (2007). Use of the weighted histogram analysis method for the analysis of simulated and parallel +tempering simulations, JCTC 3(1):26-41* + + Examples --- The example folder contains input and output files for two simple test systems: @@ -144,7 +158,6 @@ The example folder contains input and output files for two simple test systems: TODO --- -- Autocorrelation - Replica exchange License & Citing diff --git a/example/1d_cyclic/wham_uncorrelated.out b/example/1d_cyclic/wham_uncorrelated.out new file mode 100644 index 0000000..8ca4fdc --- /dev/null +++ b/example/1d_cyclic/wham_uncorrelated.out @@ -0,0 +1,101 @@ +#coord1 Free Energy +/- Probability +/- +-3.110177 7.531315 0.000000 0.003080 0.000000 +-3.047345 5.690157 0.000000 0.006443 0.000000 +-2.984513 4.243063 0.000000 0.011509 0.000000 +-2.921681 3.334686 0.000000 0.016564 0.000000 +-2.858849 2.277349 0.000000 0.025309 0.000000 +-2.796017 1.723296 0.000000 0.031604 0.000000 +-2.733186 1.246264 0.000000 0.038265 0.000000 +-2.670354 1.099867 0.000000 0.040578 0.000000 +-2.607522 0.771910 0.000000 0.046279 0.000000 +-2.544690 0.770616 0.000000 0.046303 0.000000 +-2.481858 1.265507 0.000000 0.037971 0.000000 +-2.419026 1.562335 0.000000 0.033711 0.000000 +-2.356194 1.891577 0.000000 0.029542 0.000000 +-2.293363 2.227858 0.000000 0.025816 0.000000 +-2.230531 2.488355 0.000000 0.023256 0.000000 +-2.167699 2.502265 0.000000 0.023127 0.000000 +-2.104867 2.358037 0.000000 0.024503 0.000000 +-2.042035 2.278147 0.000000 0.025301 0.000000 +-1.979203 2.974067 0.000000 0.019141 0.000000 +-1.916372 2.696600 0.000000 0.021393 0.000000 +-1.853540 2.361827 0.000000 0.024466 0.000000 +-1.790708 1.516746 0.000000 0.034332 0.000000 +-1.727876 1.526829 0.000000 0.034194 0.000000 +-1.665044 0.884114 0.000000 0.044244 0.000000 +-1.602212 0.323912 0.000000 0.055385 0.000000 +-1.539380 0.197985 0.000000 0.058253 0.000000 +-1.476549 0.000000 0.000000 0.063065 0.000000 +-1.413717 0.458247 0.000000 0.052481 0.000000 +-1.350885 1.389410 0.000000 0.036131 0.000000 +-1.288053 2.386522 0.000000 0.024225 0.000000 +-1.225221 3.743253 0.000000 0.014062 0.000000 +-1.162389 5.566654 0.000000 0.006770 0.000000 +-1.099557 7.822800 0.000000 0.002740 0.000000 +-1.036726 10.128719 0.000000 0.001087 0.000000 +-0.973894 12.199246 0.000000 0.000474 0.000000 +-0.911062 14.488129 0.000000 0.000189 0.000000 +-0.848230 16.902310 0.000000 0.000072 0.000000 +-0.785398 18.910200 0.000000 0.000032 0.000000 +-0.722566 21.241681 0.000000 0.000013 0.000000 +-0.659734 22.706373 0.000000 0.000007 0.000000 +-0.596903 24.531129 0.000000 0.000003 0.000000 +-0.534071 25.936227 0.000000 0.000002 0.000000 +-0.471239 27.000262 0.000000 0.000001 0.000000 +-0.408407 28.673293 0.000000 0.000001 0.000000 +-0.345575 29.335203 0.000000 0.000000 0.000000 +-0.282743 30.841118 0.000000 0.000000 0.000000 +-0.219911 31.983859 0.000000 0.000000 0.000000 +-0.157080 32.144015 0.000000 0.000000 0.000000 +-0.094248 33.885395 0.000000 0.000000 0.000000 +-0.031416 33.783105 0.000000 0.000000 0.000000 +0.031416 34.243727 0.000000 0.000000 0.000000 +0.094248 33.975567 0.000000 0.000000 0.000000 +0.157080 32.994787 0.000000 0.000000 0.000000 +0.219911 32.607398 0.000000 0.000000 0.000000 +0.282743 31.401902 0.000000 0.000000 0.000000 +0.345575 29.911670 0.000000 0.000000 0.000000 +0.408407 28.603574 0.000000 0.000001 0.000000 +0.471239 26.925443 0.000000 0.000001 0.000000 +0.534071 25.298070 0.000000 0.000002 0.000000 +0.596903 23.638560 0.000000 0.000005 0.000000 +0.659734 21.156231 0.000000 0.000013 0.000000 +0.722566 19.126480 0.000000 0.000029 0.000000 +0.785398 17.351953 0.000000 0.000060 0.000000 +0.848230 15.135525 0.000000 0.000146 0.000000 +0.911062 13.188112 0.000000 0.000319 0.000000 +0.973894 11.536983 0.000000 0.000618 0.000000 +1.036726 10.158328 0.000000 0.001074 0.000000 +1.099557 9.109036 0.000000 0.001636 0.000000 +1.162389 8.282343 0.000000 0.002279 0.000000 +1.225221 8.022102 0.000000 0.002530 0.000000 +1.288053 8.162415 0.000000 0.002391 0.000000 +1.350885 8.600135 0.000000 0.002006 0.000000 +1.413717 9.837348 0.000000 0.001222 0.000000 +1.476549 11.363156 0.000000 0.000663 0.000000 +1.539380 13.077849 0.000000 0.000333 0.000000 +1.602212 15.353594 0.000000 0.000134 0.000000 +1.665044 17.565051 0.000000 0.000055 0.000000 +1.727876 19.710884 0.000000 0.000023 0.000000 +1.790708 21.721260 0.000000 0.000010 0.000000 +1.853540 23.567649 0.000000 0.000005 0.000000 +1.916372 25.008817 0.000000 0.000003 0.000000 +1.979203 26.405367 0.000000 0.000002 0.000000 +2.042035 28.070821 0.000000 0.000001 0.000000 +2.104867 28.877213 0.000000 0.000001 0.000000 +2.167699 29.378146 0.000000 0.000000 0.000000 +2.230531 31.093267 0.000000 0.000000 0.000000 +2.293363 30.704994 0.000000 0.000000 0.000000 +2.356194 30.563093 0.000000 0.000000 0.000000 +2.419026 31.215952 0.000000 0.000000 0.000000 +2.481858 30.331416 0.000000 0.000000 0.000000 +2.544690 29.005123 0.000000 0.000001 0.000000 +2.607522 27.674618 0.000000 0.000001 0.000000 +2.670354 24.911788 0.000000 0.000003 0.000000 +2.733186 22.746417 0.000000 0.000007 0.000000 +2.796017 20.608288 0.000000 0.000016 0.000000 +2.858849 18.018280 0.000000 0.000046 0.000000 +2.921681 15.949215 0.000000 0.000105 0.000000 +2.984513 13.617809 0.000000 0.000268 0.000000 +3.047345 11.432193 0.000000 0.000645 0.000000 +3.110177 9.461494 0.000000 0.001420 0.000000 diff --git a/src/cli.yml b/src/cli.yml index 633adfd..7c3f714 100644 --- a/src/cli.yml +++ b/src/cli.yml @@ -101,3 +101,9 @@ args: help: Skip rows in timeseries with an index larger than this value (defaults to 1e+20) takes_value: true required: false + - uncorr: + short: g + long: uncorr + help: Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated samples (default is off). + takes_value: false + required: false diff --git a/src/correlation_analysis.rs b/src/correlation_analysis.rs new file mode 100644 index 0000000..a26bb25 --- /dev/null +++ b/src/correlation_analysis.rs @@ -0,0 +1,102 @@ +use rgsl::statistics; + +// calculates the statistical inefficiency g of the given timeseries +// the quantity g can be thought of: N/g is the number of uncorrelated +// configurations in the timeseries, where samples are separated by +// the a multiple of g +// For details, see "Chodera et al. (2007). Use of a Weighted Histogram Analysis +// Method for the Analysis of Simulated and Parallel Tempering Simulations, JCTC" +pub fn statistical_ineff(timeseries: &[f64]) -> f64 { + let n = timeseries.len(); + let autocorr = autocorrelation(timeseries); + + let mut g = 1.0; + for t in 1..(n-1) { + let c = autocorr[t-1]; + if c <= 0.0 { + break; + } + g = g + (2.0*c*(1.0-t as f64/n as f64)) + } + if g < 1.0 { + 1.0 + } else { + g + } +} + +// calculates the autocorrelation of a simeseries +fn autocorrelation(timeseries: &[f64]) -> Vec { + let n = timeseries.len(); + let mean = statistics::mean(timeseries, 1, timeseries.len()); + let d_mean = timeseries.iter().map(|x| x-mean).collect::>(); + let cov = statistics::covariance(timeseries, 1, timeseries, 1, n); + + let mut autocorr = Vec::new(); + for t in 1..(n-1) { + let tmp = d_mean[0..n-t].iter().zip(d_mean[t..n].iter()).map(|(x, y)| x*y); + let c: f64 = tmp.map(|x| x+x).sum::() / (2.0 * (n as f64-t as f64)*cov); + autocorr.push(c); + } + autocorr +} + +// The autocorrelation time of a timeseries can be deduced from the +// `statistical_ineff` by (g-1)/2.0 +pub fn autocorrelation_time(g: f64) -> f64 { + (g - 1.0) / 2.0 +} + +#[cfg(test)] +mod tests { + use std::io::{BufRead, BufReader}; + use std::fs::File; + + fn read_timeseries(filename: &str) -> Vec { + let mut timeseries: Vec = Vec::new(); + let file = File::open(filename).unwrap(); + let reader = BufReader::new(&file); + for line in reader.lines() { + let val = line.unwrap().split_whitespace() + .collect::>()[1].parse::().unwrap(); + timeseries.push(val); + } + timeseries + + } + + #[test] + fn autocorrelation() { + 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]; + for (actual, expected) in autocorr.iter().zip(expected.iter()) { + assert!((actual-expected).abs() < 0.001); + } + } + + #[test] + fn statistical_ineff() { + let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg"); + let g = super::statistical_ineff(×eries); + println!("{:?}", g); + assert!((g - 3.859).abs() < 0.001) + } + + #[test] + fn autocorrelation_time() { + let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg"); + let g = super::statistical_ineff(×eries); + let tau = super::autocorrelation_time(g); + println!("{:?}", tau); + assert!((tau - 1.430).abs() < 0.001) + } + +} \ No newline at end of file diff --git a/src/error_analysis.rs b/src/error_analysis.rs index ce8bb9e..45219e8 100644 --- a/src/error_analysis.rs +++ b/src/error_analysis.rs @@ -66,7 +66,7 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec,Ve (P_se, A_se) } -#[cfg(tests)] +#[cfg(test)] mod tests { use super::*; use super::super::k_B; @@ -99,8 +99,9 @@ mod tests { #[test] fn random_weights() { + let mut rng = StdRng::from_entropy(); let num_windows = 5; - let weights = generate_random_weights(num_windows); + let weights = generate_random_weights(num_windows, &mut rng); assert_eq!(num_windows, weights.len()); for w in weights { assert!(0.0 < w && w < 1.0); @@ -109,8 +110,9 @@ mod tests { #[test] fn random_weighted_dataset() { + let mut rng = StdRng::from_entropy(); let ds = build_hist_set(); - let rnd_weights_ds = generate_random_weighted_dataset(ds); + let rnd_weights_ds = generate_random_weighted_dataset(ds, &mut rng); println!("{:?}", rnd_weights_ds.weights); for w in rnd_weights_ds.weights { assert!(w > 0.0); diff --git a/src/io.rs b/src/io.rs index 4eba21f..6c2b904 100644 --- a/src/io.rs +++ b/src/io.rs @@ -1,6 +1,7 @@ use super::histogram::Dataset; use super::histogram::Histogram; use super::Config; +use super::correlation_analysis::{statistical_ineff, autocorrelation_time}; use std::fs::File; use std::io::prelude::*; use std::io::{BufReader,BufWriter}; @@ -115,12 +116,80 @@ 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 { - let f = File::open(window_file) +// Read a multidimensional timeseries +// The resulting vector contains one vector per dimension +fn read_timeseries(window_file: &str, cfg: &Config) -> Result>> { + let f = File::open(window_file) .chain_err(|| format!("Failed to open sample data file {}.", window_file))?; let mut buf = BufReader::new(&f); + let mut timeseries = vec![Vec::new(); cfg.dimens+1]; + + // read and parse each timeseries line + let mut line = String::new(); + let mut linecount = 0; + while buf.read_line(&mut line).chain_err(|| "Failed to read line")? > 0 { + linecount += 1; + + // skip comments and empty lines + if line.starts_with('#') || line.starts_with('@') || line.is_empty() { + line.clear(); + continue; + } + + { + let split: Vec<&str> = line.split_whitespace().collect(); + if split.len() < cfg.dimens+1 { + bail!(format!("Wrong number of columns in line {} of window file {}. Empty Line?.", linecount, window_file)); + } + + for i in 0..cfg.dimens+1 { + timeseries[i].push(split[i].parse::() + .chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))? + + ); + } + } + + line.clear(); + } + Ok(timeseries) +} + +// calculates the inefficiency for every collective variable +// filters the timeseries based on the highest inefficiency +fn uncorrelate(timeseries: Vec>, cfg: &Config) -> Vec> { + // calculate inefficiencies and find the highest one + let gs: Vec = timeseries[1..].iter().map(|ts| statistical_ineff(ts)).collect(); + let mut max_g = 1.0; + for g in gs { + if g > max_g { + max_g = g; + } + } + + // round g up + let mut trunc_g = max_g.trunc() as usize; + if (trunc_g as f64 - max_g).abs() > 0.000_000_000_1 { + trunc_g += 1; + } + + // filter correlated samples from timeseries + let prev_len = timeseries[0].len(); + let timeseries = timeseries.into_iter().map(|ts| { + ts.into_iter().step_by(trunc_g).collect::>() + }).collect::>>(); + + let new_len = timeseries[0].len(); + if cfg.verbose { + let tau = autocorrelation_time(max_g)* (timeseries[0][1]-timeseries[0][0]); + vprintln(format!("{:?}/{:?} samples are uncorrelated. {:?} samples removed from timeseries (tau={:.5})", new_len, prev_len, prev_len-new_len, tau), true); + } + timeseries +} + +// parse a time series file into a histogram +fn read_window_file(window_file: &str, cfg: &Config) -> Result { // total number of bins is the product of all dimensions length let total_bins = cfg.num_bins.iter().product(); let mut hist = vec![0.0; total_bins]; @@ -130,40 +199,26 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result { (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64) }).collect(); - // read and parse each timeseries line - let mut line = String::new(); - let mut linecount = 0; - while buf.read_line(&mut line).chain_err(|| "Failed to read line")? > 0 { - linecount += 1; - // skip comments and empty lines - if line.starts_with('#') || line.starts_with('@') || line.is_empty() { - line.clear(); - continue; + let mut timeseries: Vec> = read_timeseries(window_file, cfg)?; + + if cfg.uncorr { + timeseries = uncorrelate(timeseries, cfg); + } + + for i in 0..timeseries[0].len() { + let mut values: Vec = vec![f64::NAN; cfg.dimens+1]; + for j in 0..values.len() { + values[j] = timeseries[j][i]; } - { - - let split: Vec<&str> = line.split_whitespace().collect(); - if split.len() < cfg.dimens+1 { - bail!(format!("Wrong number of columns in line {} of window file {}. Empty Line?.", linecount, window_file)); - } - - let mut values: Vec = vec![f64::NAN; cfg.dimens+1]; - for i in 0..values.len() { - values[i] = split[i].parse::() - .chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))?; - } - - if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) { - let bin_indeces: Vec = (0..cfg.dimens).map(|dimen: usize| { - let val = values[dimen+1]; - ((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize - }).collect(); - let index = flat_index(&bin_indeces, &cfg.num_bins); - hist[index] += 1.0; - } + if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) { + let bin_indeces: Vec = (0..cfg.dimens).map(|dimen: usize| { + let val = values[dimen+1]; + ((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize + }).collect(); + let index = flat_index(&bin_indeces, &cfg.num_bins); + hist[index] += 1.0; } - line.clear(); } let num_points: f64 = hist.iter().sum(); @@ -214,7 +269,8 @@ mod tests { bootstrap: 0, bootstrap_seed: 1234, start: 0.0, - end: 1e+20 + end: 1e+20, + uncorr: false, } } @@ -233,6 +289,26 @@ mod tests { assert_approx_eq!(0.0, h.bins[7]); } + #[test] + fn read_timeseries() { + let f = "example/1d_cyclic/COLVAR+0.0.xvg"; + let cfg = cfg(); + let ts = super::read_timeseries(&f, &cfg).unwrap(); + let expected = [ + -0.153_145, + -0.377_860, + 0.010_992, + 0.123_074, + 0.108_291, + 0.261_607, + ]; + assert!(ts.len() == 2); + assert!(ts[0].len() == 5000); + println!("{:?}", ts); + for (actual, expected) in ts[1].iter().zip(expected.iter()) { + assert!((actual-expected).abs() < 0.001, format!("{:?} != {:?}", actual, expected)); + } + } #[test] fn read_data() { diff --git a/src/lib.rs b/src/lib.rs index 90ec06f..c0fbea3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -13,6 +13,7 @@ extern crate assert_approx_eq; pub mod io; pub mod histogram; pub mod error_analysis; +pub mod correlation_analysis; use histogram::Dataset; use std::f64; @@ -45,16 +46,17 @@ pub struct Config { pub bootstrap_seed: u64, pub start: f64, pub end: f64, + pub uncorr: bool, } impl fmt::Display for Config { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?}, verbose={}, tolerance={}, iterations={}, temperature={}, - cyclic={:?}, bootstrap={:?}, seed={:?}", + cyclic={:?}, uncorr={:?}, bootstrap={:?}, seed={:?}", self.metadata_file, self.hist_min, self.hist_max, self.num_bins, self.verbose, self.tolerance, self.max_iterations, self.temperature, - self.cyclic, self.bootstrap, self.bootstrap_seed) + self.cyclic, self.uncorr, self.bootstrap, self.bootstrap_seed) } } diff --git a/src/main.rs b/src/main.rs index 6764e7a..6aadf80 100644 --- a/src/main.rs +++ b/src/main.rs @@ -58,6 +58,8 @@ fn cli() -> Result { .chain_err(|| "Cannot parse start time.")?; let end: f64 = matches.value_of("end").unwrap_or("1e+20").parse() .chain_err(|| "Cannot parse end time.")?; + + let uncorr: bool = matches.is_present("uncorr"); if num_bins.len() != hist_max.len() || num_bins.len() != hist_max.len() { eprintln!("Input dimensions do not match (min: {}, max: {}, bins: {})", @@ -69,7 +71,7 @@ fn cli() -> Result { Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens, verbose, tolerance, max_iterations, temperature, cyclic, output, - bootstrap, bootstrap_seed, start, end}) + bootstrap, bootstrap_seed, start, end, uncorr}) } fn main() { diff --git a/tests/examples.rs b/tests/examples.rs index e2643f2..783b453 100644 --- a/tests/examples.rs +++ b/tests/examples.rs @@ -27,6 +27,26 @@ mod integration { assert_eq!(output_len, 0); } + #[test] + fn wham_1d_cyclic_uncorrelated() { + get_command() + .args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic", "--uncorr"]) + .args(&["--seed", "1234"]) + .args(&["-f", "example/1d_cyclic/metadata.dat"]) + .args(&["-o", "/tmp/wham_test_1d_cyclic.out"]) + .output() + .expect("failed to execute process"); + + assert!(fs::metadata("/tmp/wham_test_1d_cyclic.out").is_ok()); + let output = Command::new("diff") + .arg("/tmp/wham_test_1d_cyclic.out") + .arg("example/1d_cyclic/wham_uncorrelated.out") + .output() + .expect("failed to run diff"); + let output_len = String::from_utf8_lossy(&output.stdout).len(); + assert_eq!(output_len, 0); + } + #[test] #[ignore] // expensive fn wham_2d_cyclic() { From a856e5d34f2bffd0520c43b99e309af52ed3d73f Mon Sep 17 00:00:00 2001 From: Daniel Bauer Date: Sun, 25 Oct 2020 22:41:05 +0100 Subject: [PATCH 2/3] Update cli.yml --- src/cli.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cli.yml b/src/cli.yml index 7c3f714..15e6dad 100644 --- a/src/cli.yml +++ b/src/cli.yml @@ -1,6 +1,6 @@ name: wham version: "0.9.9" -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. From 18debf1a2862fc79b6da0f7f771f20bdb501bf95 Mon Sep 17 00:00:00 2001 From: Daniel Bauer Date: Sun, 25 Oct 2020 22:48:38 +0100 Subject: [PATCH 3/3] Update README.md --- README.md | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 3793af9..f995e79 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*).