bauer/WHAM
1
0
Fork 0
mirror of https://github.com/dnlbauer/WHAM.git synced 2026-09-11 14:45:32 +00:00
Code Issues Packages Projects Releases Wiki Activity

code formatting

Browse Source
This commit is contained in:
Daniel Bauer
2021-07-22 15:48:14 +02:00
parent ec7024a285
commit b6881a6db1
5 changed files with 347 additions and 347 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
src/error_analysis.rs
View File

@@ -73,29 +73,29 @@ mod tests {
use super::super::histogram::Histogram; use super::super::histogram::Histogram;
fn build_hist() -> Histogram { fn build_hist() -> Histogram {
Histogram::new( Histogram::new(
22, // num_points 22, // num_points
vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins
) )
} }
fn build_hist_set() -> Dataset { fn build_hist_set() -> Dataset {
let h1 = build_hist(); let h1 = build_hist();
let h2 = build_hist(); let h2 = build_hist();
let h3 = build_hist(); let h3 = build_hist();
Dataset::new( Dataset::new(
5, // num bins 5, // num bins
vec![3], vec![3],
vec![1.0], // bin width vec![1.0], // bin width
vec![0.0], // hist min vec![0.0], // hist min
vec![9.0], // hist max vec![9.0], // hist max
vec![4.5, 4.5, 4.5], // x0 vec![4.5, 4.5, 4.5], // x0
vec![10.0, 10.0, 10.0], // fc vec![10.0, 10.0, 10.0], // fc
300.0*k_B, // kT 300.0*k_B, // kT
vec![h1, h2, h3], // hists vec![h1, h2, h3], // hists
false // cyclic false // cyclic
) )
} }
#[test] #[test]
fn random_weights() { fn random_weights() {

392
src/histogram.rs
View File

@@ -3,266 +3,266 @@ use std::fmt;
// One histogram // One histogram
#[derive(Debug,Clone)] #[derive(Debug,Clone)]
pub struct Histogram { pub struct Histogram {
// total number of data points stored in the histogram // total number of data points stored in the histogram
pub num_points: u32, pub num_points: u32,
// histogram bins // histogram bins
pub bins: Vec<f64> pub bins: Vec<f64>
} }
impl Histogram { impl Histogram {
pub fn new(num_points: u32, bins: Vec<f64>) -> Histogram { pub fn new(num_points: u32, bins: Vec<f64>) -> Histogram {
Histogram {num_points, bins} Histogram {num_points, bins}
} }
} }
// a set of histograms // a set of histograms
#[derive(Debug,Clone)] #[derive(Debug,Clone)]
pub struct Dataset { pub struct Dataset {
// number of histogram windows (number of simulations) // number of histogram windows (number of simulations)
pub num_windows: usize, pub num_windows: usize,
// total number of bins // total number of bins
pub num_bins: usize, pub num_bins: usize,
// number of bins in each dimension // number of bins in each dimension
pub dimens_lengths: Vec<usize>, pub dimens_lengths: Vec<usize>,
// min values of the histogram in each dimension // min values of the histogram in each dimension
hist_min: Vec<f64>, hist_min: Vec<f64>,
// max values of the histogram in each dimension // max values of the histogram in each dimension
hist_max: Vec<f64>, hist_max: Vec<f64>,
// width of a bin in unit of its dimension // width of a bin in unit of its dimension
bin_width: Vec<f64>, bin_width: Vec<f64>,
// value of kT // value of kT
pub kT: f64, pub kT: f64,
// histogram for each window // histogram for each window
pub histograms: Vec<Histogram>, pub histograms: Vec<Histogram>,
// flag for cyclic reaction coordinates // flag for cyclic reaction coordinates
pub cyclic: bool, pub cyclic: bool,
// locations of biases // locations of biases
bias_pos: Vec<f64>, bias_pos: Vec<f64>,
// force constants of biases // force constants of biases
bias_fc: Vec<f64>, bias_fc: Vec<f64>,
// bias value cache // bias value cache
bias: Vec<f64>, bias: Vec<f64>,
// histogram weight // histogram weight
pub weights: Vec<f64>, pub weights: Vec<f64>,
} }
impl Dataset { impl Dataset {
pub fn new(num_bins: usize, dimens_lengths: Vec<usize>, bin_width: Vec<f64>, pub fn new(num_bins: usize, dimens_lengths: Vec<usize>, bin_width: Vec<f64>,
hist_min: Vec<f64>, hist_max: Vec<f64>, bias_pos: Vec<f64>, hist_min: Vec<f64>, hist_max: Vec<f64>, bias_pos: Vec<f64>,
bias_fc: Vec<f64>, kT: f64, histograms: Vec<Histogram>, cyclic: bool) -> Dataset { bias_fc: Vec<f64>, kT: f64, histograms: Vec<Histogram>, cyclic: bool) -> Dataset {
let num_windows = histograms.len(); let num_windows = histograms.len();
let bias: Vec<f64> = vec![0.0; num_bins*num_windows]; let bias: Vec<f64> = vec![0.0; num_bins*num_windows];
let weights = vec![1.0; num_windows]; let weights = vec![1.0; num_windows];
let mut ds = Dataset{ let mut ds = Dataset{
num_windows, num_windows,
num_bins, num_bins,
dimens_lengths, dimens_lengths,
bin_width, bin_width,
hist_min, hist_min,
hist_max, hist_max,
kT, kT,
histograms, histograms,
cyclic, cyclic,
bias_pos, bias_pos,
bias_fc, bias_fc,
bias, bias,
weights weights
}; };
for window in 0..num_windows { for window in 0..num_windows {
for bin in 0..num_bins { for bin in 0..num_bins {
let ndx = window * num_bins + bin; let ndx = window * num_bins + bin;
ds.bias[ndx] = ds.calc_bias(bin, window); ds.bias[ndx] = ds.calc_bias(bin, window);
} }
} }
ds ds
} }
pub fn new_weighted(ds: Dataset, weights: Vec<f64>) -> Dataset { pub fn new_weighted(ds: Dataset, weights: Vec<f64>) -> Dataset {
Dataset { Dataset {
weights, weights,
..ds ..ds
} }
} }
pub fn get_weighted_bin_count(&self, bin: usize) -> f64 { pub fn get_weighted_bin_count(&self, bin: usize) -> f64 {
self.histograms.iter().enumerate().map(|(idx,h)| self.weights[idx]*h.bins[bin]).sum() self.histograms.iter().enumerate().map(|(idx,h)| self.weights[idx]*h.bins[bin]).sum()
} }
fn expand_index(&self, bin: usize, lengths: &[usize]) -> Vec<usize> { fn expand_index(&self, bin: usize, lengths: &[usize]) -> Vec<usize> {
let mut tmp = bin; let mut tmp = bin;
let mut idx = vec![0; lengths.len()]; let mut idx = vec![0; lengths.len()];
for dimen in (1..lengths.len()).rev() { for dimen in (1..lengths.len()).rev() {
let denom: usize = lengths.iter().take(dimen).product(); let denom: usize = lengths.iter().take(dimen).product();
idx[dimen] = tmp / denom; idx[dimen] = tmp / denom;
tmp %= denom; tmp %= denom;
} }
idx[0] = tmp; idx[0] = tmp;
idx idx
} }
// get center x value for a bin // get center x value for a bin
pub fn get_coords_for_bin(&self, bin: usize) -> Vec<f64> { pub fn get_coords_for_bin(&self, bin: usize) -> Vec<f64> {
self.expand_index(bin, &self.dimens_lengths).iter().enumerate().map(|(i, dimen_bin)| { self.expand_index(bin, &self.dimens_lengths).iter().enumerate().map(|(i, dimen_bin)| {
self.hist_min[i] + self.bin_width[i]*(*dimen_bin as f64 + 0.5) self.hist_min[i] + self.bin_width[i]*(*dimen_bin as f64 + 0.5)
}).collect() }).collect()
} }
pub fn get_bias(&self, bin: usize, window: usize) -> f64 { pub fn get_bias(&self, bin: usize, window: usize) -> f64 {
let ndx = window * self.num_bins + bin; let ndx = window * self.num_bins + bin;
self.bias[ndx] self.bias[ndx]
} }
// Harmonic bias calculation: bias = 0.5*k(dx)^2 // Harmonic bias calculation: bias = 0.5*k(dx)^2
// if cyclic is true, lowest and highest bins are assumed to be // if cyclic is true, lowest and highest bins are assumed to be
// neighbors. This returns exp(U/kT) instead of U for better performance. // neighbors. This returns exp(U/kT) instead of U for better performance.
fn calc_bias(&self, bin: usize, window: usize) -> f64 { fn calc_bias(&self, bin: usize, window: usize) -> f64 {
let dimens = self.dimens_lengths.len(); let dimens = self.dimens_lengths.len();
// index of the bias value depends on the window und dimension // index of the bias value depends on the window und dimension
let bias_ndx: Vec<usize> = (0..dimens) let bias_ndx: Vec<usize> = (0..dimens)
.map(|dimen| { window * dimens + dimen }).collect(); .map(|dimen| { window * dimens + dimen }).collect();
// find the N coords, force constants and bias coords // find the N coords, force constants and bias coords
let coord = self.get_coords_for_bin(bin); let coord = self.get_coords_for_bin(bin);
let bias_fc: Vec<f64> = bias_ndx.iter().map(|ndx| { self.bias_fc[*ndx] }).collect(); let bias_fc: Vec<f64> = bias_ndx.iter().map(|ndx| { self.bias_fc[*ndx] }).collect();
let bias_pos: Vec<f64> = bias_ndx.iter().map(|ndx| { self.bias_pos[*ndx] }).collect(); let bias_pos: Vec<f64> = bias_ndx.iter().map(|ndx| { self.bias_pos[*ndx] }).collect();
let mut bias_sum = 0.0; let mut bias_sum = 0.0;
for i in 0..dimens { for i in 0..dimens {
let mut dist = (coord[i] - bias_pos[i]).abs(); let mut dist = (coord[i] - bias_pos[i]).abs();
if self.cyclic { // periodic conditions if self.cyclic { // periodic conditions
let hist_len = self.hist_max[i] - self.hist_min[i]; let hist_len = self.hist_max[i] - self.hist_min[i];
if dist > 0.5 * hist_len { if dist > 0.5 * hist_len {
dist -= hist_len; dist -= hist_len;
} }
} }
// store exp(U/kT) for better performance // store exp(U/kT) for better performance
bias_sum += 0.5 * bias_fc[i] * dist * dist bias_sum += 0.5 * bias_fc[i] * dist * dist
} }
(-bias_sum/self.kT).exp() (-bias_sum/self.kT).exp()
} }
} }
impl fmt::Display for Dataset { impl fmt::Display for Dataset {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
let mut datapoints: u32 = 0; let mut datapoints: u32 = 0;
for h in &self.histograms { for h in &self.histograms {
datapoints += h.num_points; datapoints += h.num_points;
} }
write!(f, "{} windows, {} datapoints", self.num_windows, datapoints) write!(f, "{} windows, {} datapoints", self.num_windows, datapoints)
} }
} }
#[cfg(test)] #[cfg(test)]
mod tests { mod tests {
use super::*; use super::*;
use super::super::k_B; use super::super::k_B;
macro_rules! assert_delta { macro_rules! assert_delta {
($x:expr, $y:expr, $d:expr) => { ($x:expr, $y:expr, $d:expr) => {
assert!(($x-$y).abs() < $d, "{} != {}", $x, $y) assert!(($x-$y).abs() < $d, "{} != {}", $x, $y)
} }
} }
fn build_hist() -> Histogram { fn build_hist() -> Histogram {
Histogram::new( Histogram::new(
22, // num_points 22, // num_points
vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins
) )
} }
fn build_hist_set() -> Dataset { fn build_hist_set() -> Dataset {
let h = build_hist(); let h = build_hist();
Dataset::new( Dataset::new(
5, // num bins 5, // num bins
vec![1], vec![1],
vec![1.0], // bin width vec![1.0], // bin width
vec![0.0], // hist min vec![0.0], // hist min
vec![9.0], // hist max vec![9.0], // hist max
vec![4.5], // x0 vec![4.5], // x0
vec![10.0], // fc vec![10.0], // fc
300.0*k_B, // kT 300.0*k_B, // kT
vec![h], // hists vec![h], // hists
false // cyclic false // cyclic
) )
} }
#[test] #[test]
fn calc_bias() { fn calc_bias() {
let ds = build_hist_set(); // k = 10 let ds = build_hist_set(); // k = 10
// 3th element -> x=3.5, x0=3.5 // 3th element -> x=3.5, x0=3.5
assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01); assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01);
// 8th element -> x=8.5, x0=3.5 // 8th element -> x=8.5, x0=3.5
assert_delta!(1.0, ds.calc_bias(4,0), 0.000_000_01); assert_delta!(1.0, ds.calc_bias(4,0), 0.000_000_01);
// 1st element -> x=0.5, x0=3.5. non-cyclic! // 1st element -> x=0.5, x0=3.5. non-cyclic!
assert_delta!(0.0, ds.calc_bias(0,0), 0.000_000_1); assert_delta!(0.0, ds.calc_bias(0,0), 0.000_000_1);
} }
#[test] #[test]
fn calc_biascyclic() { fn calc_biascyclic() {
let mut ds = build_hist_set(); let mut ds = build_hist_set();
ds.cyclic = true; ds.cyclic = true;
// 7th element -> x=3.5, x0=3.5 // 7th element -> x=3.5, x0=3.5
assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01); assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01);
// 8th element -> x=4.5, x0=3.5 // 8th element -> x=4.5, x0=3.5
assert_delta!(1.0, ds.calc_bias(4, 0), 0.000_000_01); assert_delta!(1.0, ds.calc_bias(4, 0), 0.000_000_01);
// 1th element -> x=0.5, x0=3.5 // 1th element -> x=0.5, x0=3.5
// cyclic flag makes bin 0 neighboring bin 9, so the distance is actually 2 // cyclic flag makes bin 0 neighboring bin 9, so the distance is actually 2
assert_delta!(0.000_000_000_000_011_776_9, ds.calc_bias(0, 0), 0.000_000_01); assert_delta!(0.000_000_000_000_011_776_9, ds.calc_bias(0, 0), 0.000_000_01);
// 2nd element -> x=1.5, x0=3.5 // 2nd element -> x=1.5, x0=3.5
assert_delta!(0.000_000_01, ds.calc_bias(1, 0), 0.000_000_01); assert_delta!(0.000_000_01, ds.calc_bias(1, 0), 0.000_000_01);
} }
#[test] #[test]
fn get_x_for_bin() { fn get_x_for_bin() {
let ds = build_hist_set(); let ds = build_hist_set();
let expected: Vec<f64> = vec![0,1,2,3,4,5,6,7,8].iter() let expected: Vec<f64> = vec![0,1,2,3,4,5,6,7,8].iter()
.map(|x| *x as f64 + 0.5).collect(); .map(|x| *x as f64 + 0.5).collect();
expected.iter().enumerate().for_each(|(i, exp)| { expected.iter().enumerate().for_each(|(i, exp)| {
assert_approx_eq!(exp, &ds.get_coords_for_bin(i)[0]); assert_approx_eq!(exp, &ds.get_coords_for_bin(i)[0]);
}) })
} }
#[test] #[test]
fn get_bin_count() { fn get_bin_count() {
let ds = Dataset::new( let ds = Dataset::new(
5, // num bins 5, // num bins
vec![1], vec![1],
vec![1.0, 1.0], // bin width vec![1.0, 1.0], // bin width
vec![0.0, 0.0], // hist min vec![0.0, 0.0], // hist min
vec![5.0, 5.0], // hist max vec![5.0, 5.0], // hist max
vec![7.5, 7.5], // x0 vec![7.5, 7.5], // x0
vec![10.0, 10.0], // fc vec![10.0, 10.0], // fc
300.0*k_B, // kT 300.0*k_B, // kT
vec![build_hist(), build_hist()], // hists vec![build_hist(), build_hist()], // hists
false // cyclic false // cyclic
); );
assert_delta!(2.0, ds.get_weighted_bin_count(0), 0.000_000_000_1); assert_delta!(2.0, ds.get_weighted_bin_count(0), 0.000_000_000_1);
assert_delta!(2.0, ds.get_weighted_bin_count(1), 0.000_000_000_1); assert_delta!(2.0, ds.get_weighted_bin_count(1), 0.000_000_000_1);
assert_delta!(6.0, ds.get_weighted_bin_count(2), 0.000_000_000_1); assert_delta!(6.0, ds.get_weighted_bin_count(2), 0.000_000_000_1);
assert_delta!(10.0, ds.get_weighted_bin_count(3), 0.000_000_000_1); assert_delta!(10.0, ds.get_weighted_bin_count(3), 0.000_000_000_1);
assert_delta!(24.0, ds.get_weighted_bin_count(4), 0.000_000_000_1); assert_delta!(24.0, ds.get_weighted_bin_count(4), 0.000_000_000_1);
} }
} }

14
src/io.rs
View File

@@ -30,8 +30,8 @@ pub fn vprintln(s: String, verbose: bool) {
// given in the metadata file. This generates at least one Dataset, // given in the metadata file. This generates at least one Dataset,
// or multiple Datasets if convdt is set in the config // or multiple Datasets if convdt is set in the config
pub fn read_data(cfg: &Config) -> Result<Vec<Dataset>> { pub fn read_data(cfg: &Config) -> Result<Vec<Dataset>> {
let mut bias_pos: Vec<f64> = Vec::new(); let mut bias_pos: Vec<f64> = Vec::new();
let mut bias_fc: Vec<f64> = Vec::new(); let mut bias_fc: Vec<f64> = Vec::new();
let mut timeseries_lengths: Vec<usize> = Vec::new(); let mut timeseries_lengths: Vec<usize> = Vec::new();
let mut paths = Vec::new(); let mut paths = Vec::new();
@@ -43,7 +43,7 @@ pub fn read_data(cfg: &Config) -> Result<Vec<Dataset>> {
// start..convdt, start..2*convdt, ... // start..convdt, start..2*convdt, ...
let mut histograms = vec![Vec::new(); dataset_boundaries.len()]; let mut histograms = vec![Vec::new(); dataset_boundaries.len()];
let kT = cfg.temperature * k_B; let kT = cfg.temperature * k_B;
let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| { let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
(cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64) (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
}).collect(); }).collect();
@@ -55,14 +55,14 @@ pub fn read_data(cfg: &Config) -> Result<Vec<Dataset>> {
// read each metadata file line and parse it // read each metadata file line and parse it
for (line_num,l) in buf.lines().enumerate() { for (line_num,l) in buf.lines().enumerate() {
let line = l.chain_err(|| "Failed to read line")?; let line = l.chain_err(|| "Failed to read line")?;
// skip comments and empty lines // skip comments and empty lines
if line.starts_with('#') || line.is_empty() { if line.starts_with('#') || line.is_empty() {
continue; continue;
} }
let split: Vec<&str> = line.split_whitespace().collect(); let split: Vec<&str> = line.split_whitespace().collect();
if split.len() < 1 + cfg.dimens * 2 { if split.len() < 1 + cfg.dimens * 2 {
bail!(format!("Wrong number of columns in line {} of metadata file. Empty Line?", line_num+1)); bail!(format!("Wrong number of columns in line {} of metadata file. Empty Line?", line_num+1));
} }

184
src/lib.rs
View File

@@ -31,18 +31,18 @@ static k_B: f64 = 0.008_314_462_1; // kJ/mol*K
// Application config // Application config
#[derive(Debug)] #[derive(Debug)]
pub struct Config { pub struct Config {
pub metadata_file: String, pub metadata_file: String,
pub hist_min: Vec<f64>, pub hist_min: Vec<f64>,
pub hist_max: Vec<f64>, pub hist_max: Vec<f64>,
pub num_bins: Vec<usize>, pub num_bins: Vec<usize>,
pub dimens: usize, pub dimens: usize,
pub verbose: bool, pub verbose: bool,
pub tolerance: f64, pub tolerance: f64,
pub max_iterations: usize, pub max_iterations: usize,
pub temperature: f64, pub temperature: f64,
pub cyclic: bool, pub cyclic: bool,
pub output: String, pub output: String,
pub bootstrap: usize, pub bootstrap: usize,
pub bootstrap_seed: u64, pub bootstrap_seed: u64,
pub start: f64, pub start: f64,
pub end: f64, pub end: f64,
@@ -52,7 +52,7 @@ pub struct Config {
} }
impl fmt::Display for Config { impl fmt::Display for Config {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?}, write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?},
verbose={}, tolerance={}, iterations={}, temperature={}, verbose={}, tolerance={}, iterations={}, temperature={},
cyclic={:?}, uncorr={:?}, bootstrap={:?}, seed={:?}, cyclic={:?}, uncorr={:?}, bootstrap={:?}, seed={:?},
@@ -70,7 +70,7 @@ impl fmt::Display for Config {
fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool { fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool {
// calculates abs diff between every old and new F and checks if any // calculates abs diff between every old and new F and checks if any
// is larger than tolerance // is larger than tolerance
!new_F.iter() !new_F.iter()
.zip(old_F.iter()) .zip(old_F.iter())
.map(|x| { (x.0-x.1).abs() }) .map(|x| { (x.0-x.1).abs() })
.any(|diff| { diff > tolerance }) .any(|diff| { diff > tolerance })
@@ -81,13 +81,13 @@ fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool {
// P(x) = \frac {\sum_{i=1}^N{n_i(x)}} // P(x) = \frac {\sum_{i=1}^N{n_i(x)}}
// {\sum_{i=1}^N{ N_i exp(\beta [F_i - U_{bias,i}(x)])}} // {\sum_{i=1}^N{ N_i exp(\beta [F_i - U_{bias,i}(x)])}}
fn calc_bin_probability(bin: usize, dataset: &Dataset, F: &[f64]) -> f64 { fn calc_bin_probability(bin: usize, dataset: &Dataset, F: &[f64]) -> f64 {
let mut denom_sum: f64 = 0.0; let mut denom_sum: f64 = 0.0;
let bin_count: f64 = dataset.get_weighted_bin_count(bin); let bin_count: f64 = dataset.get_weighted_bin_count(bin);
for (window, h) in dataset.histograms.iter().enumerate() { for (window, h) in dataset.histograms.iter().enumerate() {
let bias = dataset.get_bias(bin, window); let bias = dataset.get_bias(bin, window);
denom_sum += (dataset.weights[window] * h.num_points as f64) denom_sum += (dataset.weights[window] * h.num_points as f64)
* bias * F[window]; * bias * F[window];
} }
bin_count / denom_sum bin_count / denom_sum
} }
@@ -109,26 +109,26 @@ fn calc_window_F(window: usize, dataset: &Dataset, P: &[f64]) -> f64 {
// offsets F based on previous bias offsets F_prev. This updates the values in // offsets F based on previous bias offsets F_prev. This updates the values in
// vectors F and P. // vectors F and P.
fn perform_wham_iteration(dataset: &Dataset, F_prev: &[f64], F: &mut Vec<f64>, P: &mut Vec<f64>) { fn perform_wham_iteration(dataset: &Dataset, F_prev: &[f64], F: &mut Vec<f64>, P: &mut Vec<f64>) {
// Update P // Update P
// evaluate first WHAM equation for each bin to // evaluate first WHAM equation for each bin to
// estimate probabilities based on previous offsets (F_prev)) // estimate probabilities based on previous offsets (F_prev))
(0..dataset.num_bins).into_par_iter() (0..dataset.num_bins).into_par_iter()
.map(|bin| { calc_bin_probability(bin, dataset, F_prev) }) .map(|bin| { calc_bin_probability(bin, dataset, F_prev) })
.collect_into_vec(P); .collect_into_vec(P);
// Update F // Update F
// evaluate second WHAM equation for each window to // evaluate second WHAM equation for each window to
// estimate new bias offsets from propabilities // estimate new bias offsets from propabilities
(0..dataset.num_windows).into_par_iter() (0..dataset.num_windows).into_par_iter()
.map(|window| {calc_window_F(window, dataset, P)} ) .map(|window| {calc_window_F(window, dataset, P)} )
.collect_into_vec(F); .collect_into_vec(F);
} }
// Full WHAM calculation. Calls `perform_wham_iteration` until convergence // Full WHAM calculation. Calls `perform_wham_iteration` until convergence
// criteria are met or max iterations reached. // criteria are met or max iterations reached.
pub fn perform_wham(cfg: &Config, dataset: &Dataset) pub fn perform_wham(cfg: &Config, dataset: &Dataset)
-> Result<(Vec<f64>, Vec<f64>, Vec<f64>)> { -> Result<(Vec<f64>, Vec<f64>, Vec<f64>)> {
// allocate required vectors. // allocate required vectors.
// bin probability // bin probability
let mut P: Vec<f64> = vec![f64::NAN; dataset.num_bins]; let mut P: Vec<f64> = vec![f64::NAN; dataset.num_bins];
@@ -177,10 +177,10 @@ pub fn perform_wham(cfg: &Config, dataset: &Dataset)
} }
if iteration == cfg.max_iterations { if iteration == cfg.max_iterations {
bail!("WHAM not converged! (max iterations reached)"); bail!("WHAM not converged! (max iterations reached)");
} }
Ok((P, F, F_prev)) Ok((P, F, F_prev))
} }
pub fn run(cfg: &Config) -> Result<()>{ pub fn run(cfg: &Config) -> Result<()>{
@@ -227,17 +227,17 @@ pub fn run(cfg: &Config) -> Result<()>{
// get average difference between two bias offset sets // get average difference between two bias offset sets
fn diff_avg(F: &[f64], F_prev: &[f64]) -> f64 { fn diff_avg(F: &[f64], F_prev: &[f64]) -> f64 {
let mut F_sum: f64 = 0.0; let mut F_sum: f64 = 0.0;
for i in 0..F.len() { for i in 0..F.len() {
F_sum += (F[i]-F_prev[i]).abs() F_sum += (F[i]-F_prev[i]).abs()
} }
F_sum / F.len() as f64 F_sum / F.len() as f64
} }
// calculate the normalized free energy from probability values // calculate the normalized free energy from probability values
fn calc_free_energy(dataset: &Dataset, P: &[f64]) -> Vec<f64> { fn calc_free_energy(dataset: &Dataset, P: &[f64]) -> Vec<f64> {
let mut minimum = f64::MAX; let mut minimum = f64::MAX;
let mut free_energy: Vec<f64> = P.iter() let mut free_energy: Vec<f64> = P.iter()
.map(|p| { .map(|p| {
-dataset.kT * p.ln() -dataset.kT * p.ln()
}) })
@@ -257,28 +257,28 @@ fn calc_free_energy(dataset: &Dataset, P: &[f64]) -> Vec<f64> {
// Print the current WHAM iteration state. Dumps the PMF and associated vectors // Print the current WHAM iteration state. Dumps the PMF and associated vectors
fn dump_state(dataset: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], fn dump_state(dataset: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64],
P_std: &[f64], A: &[f64], A_std: &[f64]) { P_std: &[f64], A: &[f64], A_std: &[f64]) {
// TODO fix output of F/F_prev // TODO fix output of F/F_prev
let out = std::io::stdout(); let out = std::io::stdout();
let mut lock = out.lock(); let mut lock = out.lock();
writeln!(lock, "# PMF").unwrap(); writeln!(lock, "# PMF").unwrap();
writeln!(lock, "#bin\t\tFree Energy\t\t+/-\t\tP(x)\t\t+/-").unwrap(); writeln!(lock, "#bin\t\tFree Energy\t\t+/-\t\tP(x)\t\t+/-").unwrap();
for bin in 0..dataset.num_bins { for bin in 0..dataset.num_bins {
writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}", writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}",
bin, A[bin], A_std[bin], P[bin], P_std[bin]).unwrap(); bin, A[bin], A_std[bin], P[bin], P_std[bin]).unwrap();
} }
writeln!(lock, "# Bias offsets").unwrap(); writeln!(lock, "# Bias offsets").unwrap();
writeln!(lock, "#Window\t\tF\t\tF_prev").unwrap(); writeln!(lock, "#Window\t\tF\t\tF_prev").unwrap();
for window in 0..dataset.num_windows { for window in 0..dataset.num_windows {
writeln!(lock, "{}\t{:9.5}\t{:8.8}", writeln!(lock, "{}\t{:9.5}\t{:8.8}",
window, F[window], (F[window]-F_prev[window]).abs()).unwrap(); window, F[window], (F[window]-F_prev[window]).abs()).unwrap();
} }
} }
#[cfg(test)] #[cfg(test)]
mod tests { mod tests {
use super::histogram::{Dataset,Histogram}; use super::histogram::{Dataset,Histogram};
use std::f64; use std::f64;
use super::k_B; use super::k_B;
macro_rules! assert_delta { macro_rules! assert_delta {
@@ -288,65 +288,65 @@ mod tests {
} }
fn create_test_dataset() -> Dataset { fn create_test_dataset() -> Dataset {
let h1 = Histogram::new(10, vec![0.0, 1.0, 1.0, 8.0, 0.0]); let h1 = Histogram::new(10, vec![0.0, 1.0, 1.0, 8.0, 0.0]);
let h2 = Histogram::new(10, vec![0.0, 0.0, 8.0, 1.0, 1.0]); let h2 = Histogram::new(10, vec![0.0, 0.0, 8.0, 1.0, 1.0]);
Dataset::new(5, vec![5], vec![1.0], vec![0.0], vec![4.0], Dataset::new(5, vec![5], vec![1.0], vec![0.0], vec![4.0],
vec![1.0, 1.0], vec![10.0, 10.0], 300.0*k_B, vec![h1, h2], false) vec![1.0, 1.0], vec![10.0, 10.0], 300.0*k_B, vec![h1, h2], false)
} }
#[test]
fn is_converged() {
let new = vec![1.0,1.0];
let old = vec![0.95, 1.0];
let tolerance = 0.1;
let converged = super::is_converged(&old, &new, tolerance);
assert!(converged);
let old = vec![0.8, 1.0];
let converged = super::is_converged(&old, &new, tolerance);
assert!(!converged);
}
#[test] #[test]
fn calc_bin_probability() { fn is_converged() {
let dataset = create_test_dataset(); let new = vec![1.0,1.0];
let F = vec![1.0; dataset.num_bins] ; let old = vec![0.95, 1.0];
let tolerance = 0.1;
let converged = super::is_converged(&old, &new, tolerance);
assert!(converged);
let old = vec![0.8, 1.0];
let converged = super::is_converged(&old, &new, tolerance);
assert!(!converged);
}
#[test]
fn calc_bin_probability() {
let dataset = create_test_dataset();
let F = vec![1.0; dataset.num_bins] ;
let expected = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93, let expected = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93,
124_226.700_033_77, 2_308_526_035.528_374_7); 124_226.700_033_77, 2_308_526_035.528_374_7);
expected.iter().enumerate().for_each(|(i, exp)| { expected.iter().enumerate().for_each(|(i, exp)| {
let p = super::calc_bin_probability(i, &dataset, &F); let p = super::calc_bin_probability(i, &dataset, &F);
assert_delta!(exp, p, 0.000_000_1); assert_delta!(exp, p, 0.000_000_1);
}) })
} }
#[test] #[test]
fn calc_bias_offset() { fn calc_bias_offset() {
let dataset = create_test_dataset(); let dataset = create_test_dataset();
let probability = vec!(0.0, 0.1, 0.2, 0.3, 0.4); let probability = vec!(0.0, 0.1, 0.2, 0.3, 0.4);
let expected = vec!(15.927_477_169_990_633, 15.927_477_169_990_633); let expected = vec!(15.927_477_169_990_633, 15.927_477_169_990_633);
expected.iter().enumerate().for_each(|(i, exp)| { expected.iter().enumerate().for_each(|(i, exp)| {
let F = super::calc_window_F(i, &dataset, &probability); let F = super::calc_window_F(i, &dataset, &probability);
assert_delta!(exp, F, 0.000_000_1); assert_delta!(exp, F, 0.000_000_1);
}) })
} }
#[test] #[test]
fn perform_wham_iteration() { fn perform_wham_iteration() {
let dataset = create_test_dataset(); let dataset = create_test_dataset();
let prev_F = vec![1.0; dataset.num_windows]; let prev_F = vec![1.0; dataset.num_windows];
let mut F = vec![f64::NAN; dataset.num_windows]; let mut F = vec![f64::NAN; dataset.num_windows];
let mut P = vec![f64::NAN; dataset.num_bins]; let mut P = vec![f64::NAN; dataset.num_bins];
super::perform_wham_iteration(&dataset, &prev_F, &mut F, &mut P); super::perform_wham_iteration(&dataset, &prev_F, &mut F, &mut P);
let expected_F = vec!(1.0, 1.0); let expected_F = vec!(1.0, 1.0);
let expected_P = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93, let expected_P = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93,
124_226.700_033_77, 2_308_526_035.528_374_7); 124_226.700_033_77, 2_308_526_035.528_374_7);
for bin in 0..dataset.num_bins { for bin in 0..dataset.num_bins {
assert_delta!(expected_P[bin], P[bin], 0.01) assert_delta!(expected_P[bin], P[bin], 0.01)
} }
for window in 0..dataset.num_windows { for window in 0..dataset.num_windows {
assert_delta!(expected_F[window], F[window], 0.01) assert_delta!(expected_F[window], F[window], 0.01)
} }
} }
} }

54
src/main.rs
View File

@@ -13,20 +13,20 @@ use std::process;
// Parse command line arguments into a Config struct // Parse command line arguments into a Config struct
fn cli() -> Result<Config> { fn cli() -> Result<Config> {
let yaml = load_yaml!("cli.yml"); let yaml = load_yaml!("cli.yml");
let matches = App::from_yaml(yaml).get_matches(); let matches = App::from_yaml(yaml).get_matches();
let metadata_file = matches.value_of("metadata").unwrap().to_string(); let metadata_file = matches.value_of("metadata").unwrap().to_string();
let verbose: bool = matches.is_present("verbose"); let verbose: bool = matches.is_present("verbose");
let temperature: f64 = matches.value_of("temperature").unwrap().parse() let temperature: f64 = matches.value_of("temperature").unwrap().parse()
.chain_err(|| "Cannot read temperature.")?; .chain_err(|| "Cannot read temperature.")?;
let tolerance: f64 = matches.value_of("tolerance").unwrap_or("0.000001").parse() let tolerance: f64 = matches.value_of("tolerance").unwrap_or("0.000001").parse()
.chain_err(|| "Cannot read tolerance.")?; .chain_err(|| "Cannot read tolerance.")?;
let max_iterations: usize = matches.value_of("iterations").unwrap_or("100000").parse() let max_iterations: usize = matches.value_of("iterations").unwrap_or("100000").parse()
.chain_err(|| "Cannot parse iterations.")?; .chain_err(|| "Cannot parse iterations.")?;
let output = matches.value_of("output").unwrap_or("wham.out").to_string(); let output = matches.value_of("output").unwrap_or("wham.out").to_string();
let cyclic: bool = matches.is_present("cyclic"); let cyclic: bool = matches.is_present("cyclic");
let hist_min: Vec<f64> = matches.value_of("min_hist").unwrap() let hist_min: Vec<f64> = matches.value_of("min_hist").unwrap()
.split(',').map(|x| { .split(',').map(|x| {
if x.to_ascii_lowercase() == "pi" { if x.to_ascii_lowercase() == "pi" {
std::f64::consts::PI std::f64::consts::PI
@@ -36,7 +36,7 @@ fn cli() -> Result<Config> {
x.parse().unwrap() x.parse().unwrap()
} }
}).collect(); }).collect();
let hist_max: Vec<f64> = matches.value_of("max_hist").unwrap() let hist_max: Vec<f64> = matches.value_of("max_hist").unwrap()
.split(',').map(|x| { .split(',').map(|x| {
if x.to_ascii_lowercase() == "pi" { if x.to_ascii_lowercase() == "pi" {
std::f64::consts::PI std::f64::consts::PI
@@ -46,10 +46,10 @@ fn cli() -> Result<Config> {
x.parse().unwrap() x.parse().unwrap()
} }
}).collect(); }).collect();
let num_bins: Vec<usize> = matches.value_of("bins").unwrap() let num_bins: Vec<usize> = matches.value_of("bins").unwrap()
.split(',').map(|x| { x.parse().unwrap() }).collect(); .split(',').map(|x| { x.parse().unwrap() }).collect();
let bootstrap: usize = matches.value_of("bootstrap").unwrap_or("0").parse() let bootstrap: usize = matches.value_of("bootstrap").unwrap_or("0").parse()
.chain_err(|| "Cannot parse bootstrap iteration.")?; .chain_err(|| "Cannot parse bootstrap iteration.")?;
let bootstrap_seed: u64 = matches.value_of("bootstrap_seed") let bootstrap_seed: u64 = matches.value_of("bootstrap_seed")
.unwrap_or({ .unwrap_or({
let mut rng = rand::thread_rng(); let mut rng = rand::thread_rng();
@@ -78,20 +78,20 @@ fn cli() -> Result<Config> {
let ignore_empty: bool = matches.is_present("ignore_empty"); let ignore_empty: bool = matches.is_present("ignore_empty");
Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens, Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens,
verbose, tolerance, max_iterations, temperature, cyclic, output, verbose, tolerance, max_iterations, temperature, cyclic, output,
bootstrap, bootstrap_seed, start, end, uncorr, convdt, ignore_empty}) bootstrap, bootstrap_seed, start, end, uncorr, convdt, ignore_empty})
} }
fn main() { fn main() {
let cfg = cli().expect("Failed to parse CLI."); let cfg = cli().expect("Failed to parse CLI.");
if let Err(error) = wham::run(&cfg) { if let Err(error) = wham::run(&cfg) {
eprintln!("Error: {}", error); eprintln!("Error: {}", error);
for e in error.iter().skip(1) { for e in error.iter().skip(1) {
eprintln!("Reason: {}", e) eprintln!("Reason: {}", e)
} }
process::exit(1); process::exit(1);
} }
} }