diff --git a/src/histogram.rs b/src/histogram.rs index 2be9dac..9504726 100644 --- a/src/histogram.rs +++ b/src/histogram.rs @@ -14,18 +14,18 @@ pub struct Histogram { pub num_points: u32, // histogram bins - pub bins: Vec + pub bins: Vec } impl Histogram { - pub fn new(first: usize, last: usize, num_points: u32, bins: Vec) -> Histogram { + pub fn new(first: usize, last: usize, num_points: u32, bins: Vec) -> Histogram { assert_eq!(last-first+1, bins.len(), "histogram length does not match first/last."); Histogram {first, last, num_points, bins} } // Returns the value of a bin if the bin is present in this // histogram - pub fn get_bin_count(&self, bin: usize) -> Option { + pub fn get_bin_count(&self, bin: usize) -> Option { if bin < self.first || bin > self.last { None } else { @@ -44,16 +44,16 @@ pub struct Dataset { pub num_bins: usize, // min value of the histogram - pub hist_min: f32, + pub hist_min: f64, // max value of the histogram - pub hist_max: f32, + pub hist_max: f64, // width of a bin in unit of x - pub bin_width: f32, + pub bin_width: f64, // value of kT - pub kT: f32, + pub kT: f64, // histogram for each window pub histograms: Vec, @@ -62,20 +62,20 @@ pub struct Dataset { pub cyclic: bool, // locations of biases - bias_x0: Vec, + bias_x0: Vec, // force constants of biases - bias_fc: Vec, + bias_fc: Vec, // bias value cache - bias: RefCell>> + bias: RefCell>> } impl Dataset { - pub fn new(num_bins: usize, bin_width: f32, hist_min: f32, hist_max: f32, bias_x0: Vec, bias_fc: Vec, kT: f32, histograms: Vec, cyclic: bool) -> Dataset { + pub fn new(num_bins: usize, bin_width: f64, hist_min: f64, hist_max: f64, bias_x0: Vec, bias_fc: Vec, kT: f64, histograms: Vec, cyclic: bool) -> Dataset { let num_windows = histograms.len(); - let bias: RefCell>> = RefCell::new(vec![None; num_bins*num_windows]); + let bias: RefCell>> = RefCell::new(vec![None; num_bins*num_windows]); Dataset{ num_windows, num_bins, @@ -95,7 +95,7 @@ impl Dataset { // Harmonic bias calculation: bias = 0.5*k(dx)^2 // if cyclic is true, lowest and highest bins are assumed to be // neighbors - pub fn calc_bias(&self, bin: usize, window: usize) -> f32 { + pub fn calc_bias(&self, bin: usize, window: usize) -> f64 { let ndx = bin + (self.num_bins*window); let mut cache = self.bias.borrow_mut(); match cache[ndx] { @@ -117,8 +117,8 @@ impl Dataset { } // get center x value for a bin - pub fn get_x_for_bin(&self, bin: usize) -> f32 { - self.hist_min + self.bin_width * ((bin as f32) + 0.5) + pub fn get_x_for_bin(&self, bin: usize) -> f64 { + self.hist_min + self.bin_width * ((bin as f64) + 0.5) } } @@ -212,8 +212,8 @@ mod tests { #[test] fn get_x_for_bin() { let ds = build_hist_set(); - let expected: Vec = vec![0,1,2,3,4,5,6,7,8].iter() - .map(|x| *x as f32 + 0.5).collect(); + let expected: Vec = vec![0,1,2,3,4,5,6,7,8].iter() + .map(|x| *x as f64 + 0.5).collect(); for i in 0..9 { assert_eq!(expected[i], ds.get_x_for_bin(i)); } diff --git a/src/io.rs b/src/io.rs index da39a45..35aae95 100644 --- a/src/io.rs +++ b/src/io.rs @@ -29,8 +29,8 @@ pub fn vprintln(s: String, verbose: bool) { // Read input data into a histogram set by iterating over input files // given in the metadata file pub fn read_data(cfg: &Config) -> Option { - let mut bias_x0: Vec = Vec::new(); - let mut bias_fc: Vec = Vec::new(); + let mut bias_x0: Vec = Vec::new(); + let mut bias_fc: Vec = Vec::new(); let mut histograms: Vec = Vec::new(); let kT = cfg.temperature * k_B; let f = File::open(&cfg.metadata_file).unwrap_or_else(|x| { @@ -45,7 +45,7 @@ pub fn read_data(cfg: &Config) -> Option { if line.starts_with("#") || line.len() == 0 { continue; } - let (path, x0, k) = scan_fmt!(&line, "{} {} {}", String, f32, f32); + let (path, x0, k) = scan_fmt!(&line, "{} {} {}", String, f64, f64); let path = get_relative_path(&cfg.metadata_file, &path.unwrap()); match read_window_file(&path, cfg) { @@ -72,7 +72,7 @@ pub fn read_data(cfg: &Config) -> Option { } if histograms.len() > 0 { - let bin_width = (cfg.hist_max - cfg.hist_min)/(cfg.num_bins as f32); + let bin_width = (cfg.hist_max - cfg.hist_min)/(cfg.num_bins as f64); Some(Dataset::new(cfg.num_bins, bin_width, cfg.hist_min, cfg.hist_max, bias_x0, bias_fc, kT, histograms, cfg.cyclic)) } else { None @@ -88,7 +88,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option { let buf = BufReader::new(&f); let mut global_hist = vec![0.0; cfg.num_bins]; - let bin_width = (cfg.hist_max - cfg.hist_min)/(cfg.num_bins as f32); + let bin_width = (cfg.hist_max - cfg.hist_min)/(cfg.num_bins as f64); for l in buf.lines() { let line = l.unwrap(); // skip comments and empty lines @@ -96,7 +96,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option { continue; } - let (_, x) = scan_fmt!(&line, "{} {}", f32, f32); + let (_, x) = scan_fmt!(&line, "{} {}", f64, f64); match x { Some(x) => { @@ -130,12 +130,12 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option { global_hist.truncate(max_bin+1); global_hist.drain(..min_bin); - let num_points: f32 = global_hist.iter().sum(); + let num_points: f64 = global_hist.iter().sum(); Some(Histogram::new(min_bin, max_bin, num_points as u32, global_hist)) } } -pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec, prob: &Vec) -> Result<(), Box> { +pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec, prob: &Vec) -> Result<(), Box> { let mut output = File::create(out_file)?; writeln!(output, "#{:8}\t{:8}\t{:8}", "x", "Free Energy", "Probability"); for bin in 0..free.len() { @@ -191,7 +191,7 @@ mod tests { assert_eq!(cfg.num_bins, ds.num_bins); assert_eq!(cfg.hist_min, ds.hist_min); assert_eq!(cfg.hist_max, ds.hist_max); - let expected_bin_width = (cfg.hist_max - cfg.hist_min)/cfg.num_bins as f32; + let expected_bin_width = (cfg.hist_max - cfg.hist_min)/cfg.num_bins as f64; assert_eq!(expected_bin_width, ds.bin_width); // bias fields are private // assert_eq!(vec![0.0, 1.0], ds.bias_x0); diff --git a/src/lib.rs b/src/lib.rs index a089daf..ebc7b61 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -9,23 +9,23 @@ pub mod histogram; use std::error::Error; use std::result::Result; use histogram::{Dataset,Histogram}; -use std::f32; +use std::f64; use std::fmt; #[allow(non_upper_case_globals)] -static k_B: f32 = 0.0083144621; // kJ/mol*K +static k_B: f64 = 0.0083144621; // kJ/mol*K // Application config #[derive(Debug)] pub struct Config { pub metadata_file: String, - pub hist_min: f32, - pub hist_max: f32, + pub hist_min: f64, + pub hist_max: f64, pub num_bins: usize, pub verbose: bool, - pub tolerance: f32, + pub tolerance: f64, pub max_iterations: usize, - pub temperature: f32, + pub temperature: f64, pub cyclic: bool, pub output: String, } @@ -41,7 +41,7 @@ impl fmt::Display for Config { // Checks for convergence between two WHAM iterations. WHAM is considered as // converged if the absolute difference for the calculated bias offset is // smaller then a tolerance value for every simulation window. -fn is_converged(old_F: &Vec, new_F: &Vec, tolerance: f32) -> bool { +fn is_converged(old_F: &Vec, new_F: &Vec, tolerance: f64) -> bool { !new_F.iter().zip(old_F.iter()) .map(|x| { (x.0-x.1).abs() }) .any(|diff| { diff > tolerance }) @@ -49,9 +49,9 @@ fn is_converged(old_F: &Vec, new_F: &Vec, tolerance: f32) -> bool { // estimate the probability of a bin of the histogram set based on F values // This evaluates the first WHAM equation for each bin -fn calc_bin_probability(bin: usize, ds: &Dataset, F: &Vec) -> f32 { - let mut denom_sum = 0.0; - let mut bin_count = 0.0; +fn calc_bin_probability(bin: usize, ds: &Dataset, F: &Vec) -> f64 { + let mut denom_sum: f64 = 0.0; + let mut bin_count: f64 = 0.0; for window in 0..ds.num_windows { let h: &Histogram = &ds.histograms[window]; if let Some(count) = h.get_bin_count(bin) { @@ -59,16 +59,16 @@ fn calc_bin_probability(bin: usize, ds: &Dataset, F: &Vec) -> f32 { } let bias = ds.calc_bias(bin, window); let bias_offset = ((F[window] - bias) / ds.kT).exp(); - denom_sum += (h.num_points as f32) * bias_offset; + denom_sum += (h.num_points as f64) * bias_offset; } bin_count / denom_sum } // estimate the bias offset F of the histogram based on given probabilities // This evaluates the second WHAM equation for each window -fn calc_window_F(window: usize, ds: &Dataset, P: &Vec) -> f32 { - let bf_sum: f32 = (0..ds.num_bins).zip(P.iter()) // zip bins and P - .map(|x: (usize, &f32)| { +fn calc_window_F(window: usize, ds: &Dataset, P: &Vec) -> f64 { + let bf_sum: f64 = (0..ds.num_bins).zip(P.iter()) // zip bins and P + .map(|x: (usize, &f64)| { x.1 * (-ds.calc_bias(x.0, window)/ds.kT).exp() }).sum(); -ds.kT * bf_sum.ln() @@ -77,7 +77,7 @@ fn calc_window_F(window: usize, ds: &Dataset, P: &Vec) -> f32 { // One full WHAM iteration includes calculation of new probabilities P and // new bias offsets F based on previous bias offsets F_prev. This updates // the values in vectors F and P -fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec,F: &mut Vec, P: &mut Vec) { +fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec,F: &mut Vec, P: &mut Vec) { // reset bias offsets for window in 0..ds.num_windows { F[window] = 0.0; @@ -98,7 +98,7 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec,F: &mut Vec, P: & // ds.bias_x0[window], // x); // let bf = ((F_prev[window]-bias) / ds.kT).exp(); - // denom += ds.histograms[window].num_points as f32* bf + // denom += ds.histograms[window].num_points as f64* bf // } // P[bin] = num / denom; @@ -141,18 +141,18 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec,F: &mut Vec, P: & } // get average difference between two bias offset sets -fn diff_avg(F: &Vec, F_prev: &Vec) -> f32 { - let mut F_sum = 0.0; +fn diff_avg(F: &Vec, F_prev: &Vec) -> f64 { + let mut F_sum: f64 = 0.0; for i in 0..F.len() { F_sum += (F[i]-F_prev[i]).abs() } - F_sum / F.len() as f32 + F_sum / F.len() as f64 } // calculate the normalized free energy from normalized probability values -fn free_energy(ds: &Dataset, P: &mut Vec, A: &mut Vec) { - let mut bin_min = f32::MAX; +fn free_energy(ds: &Dataset, P: &mut Vec, A: &mut Vec) { + let mut bin_min = f64::MAX; // Free energy calculation for bin in 0..ds.num_bins { @@ -190,10 +190,10 @@ pub fn run(cfg: &Config) -> Result<(), Box>{ println!("{}",&histograms); // allocate only once for better performance - let mut F_prev = vec![f32::INFINITY; histograms.num_windows]; - let mut F = vec![0.0; histograms.num_windows]; - let mut P = vec![f32::NAN; histograms.num_bins]; - let mut A = vec![f32::NAN; histograms.num_bins]; + let mut F_prev: Vec = vec![f64::INFINITY; histograms.num_windows]; + let mut F: Vec = vec![0.0; histograms.num_windows]; + let mut P: Vec = vec![f64::NAN; histograms.num_bins]; + let mut A: Vec = vec![f64::NAN; histograms.num_bins]; // perform WHAM until convergence let mut iteration = 0; @@ -240,7 +240,7 @@ pub fn run(cfg: &Config) -> Result<(), Box>{ Ok(()) } -fn dump_state(ds: &Dataset, F: &Vec, F_prev: &Vec, P: &Vec, A: &Vec) { +fn dump_state(ds: &Dataset, F: &Vec, F_prev: &Vec, P: &Vec, A: &Vec) { println!("# PMF"); println!("#x\t\tFree Energy\t\tP(x)"); for bin in 0..ds.num_bins { @@ -257,7 +257,7 @@ fn dump_state(ds: &Dataset, F: &Vec, F_prev: &Vec, P: &Vec, A: &V #[cfg(test)] mod tests { use super::histogram::{Dataset,Histogram}; - use std::f32; + use std::f64; #[test] fn is_converged() { @@ -278,7 +278,7 @@ mod tests { Dataset::new(4, 1.0, 0.0, 4.0, vec![1.0, 2.0], vec![10.0, 10.0], 2.479, vec![h1, h2], false) } - fn assert_near(a: f32, b: f32, tolerance: f32) { + fn assert_near(a: f64, b: f64, tolerance: f64) { let d = (a-b).abs(); assert!(d <= tolerance, "Values are not close: {}, {}, d={}", &a, &b, &d); } @@ -317,9 +317,9 @@ mod tests { let ds = create_test_ds(); let prev_F = vec![0.0; ds.num_windows]; let mut F = vec![0.0; ds.num_windows]; - let mut P = vec![f32::NAN; ds.num_bins]; + let mut P = vec![f64::NAN; ds.num_bins]; super::perform_wham_iteration(&ds, &prev_F, &mut F, &mut P); - let expected_F = vec!(0.0, -0.846); + let expected_F = vec!(0.5948, -0.2513); let expected_P = vec!(0.959, 0.331, 0.656, 46.750); for bin in 0..ds.num_bins { assert_near(expected_P[bin], P[bin], 0.01) diff --git a/src/main.rs b/src/main.rs index 7c486ce..c4b1844 100644 --- a/src/main.rs +++ b/src/main.rs @@ -13,12 +13,12 @@ fn cli() -> Result> { let yaml = load_yaml!("cli.yml"); let matches = App::from_yaml(yaml).get_matches(); let metadata_file = matches.value_of("metadata").unwrap().to_string(); - let hist_min: f32 = matches.value_of("min_hist").unwrap().parse()?; - let hist_max: f32 = matches.value_of("max_hist").unwrap().parse()?; + let hist_min: f64 = matches.value_of("min_hist").unwrap().parse()?; + let hist_max: f64 = matches.value_of("max_hist").unwrap().parse()?; let num_bins: usize = matches.value_of("bins").unwrap().parse()?; let verbose: bool = matches.is_present("verbose"); - let temperature: f32 = matches.value_of("temperature").unwrap().parse()?; - let tolerance: f32 = matches.value_of("tolerance").unwrap_or("0.000001").parse()?; + let temperature: f64 = matches.value_of("temperature").unwrap().parse()?; + let tolerance: f64 = matches.value_of("tolerance").unwrap_or("0.000001").parse()?; let max_iterations: usize = matches.value_of("iterations").unwrap_or("100000").parse()?; let cyclic: bool = matches.is_present("cyclic"); let output = matches.value_of("output").unwrap_or("wham.out").to_string(); diff --git a/tests/data/window_0.0.dat b/tests/data/window_0.0.dat index 40cb7b0..483f1cf 100644 --- a/tests/data/window_0.0.dat +++ b/tests/data/window_0.0.dat @@ -1,12 +1,12 @@ -1 0.1 +1 0.100000001 2 0.15 -3 0.2 +3 0.200000001 4 0.25 -5 0.3 +5 0.300000001 6 0.35 -7 0.4 +7 0.400000001 8 0.45 -9 0.5 +9 0.500000001 1 0.55 -11 0.6 -12 99 \ No newline at end of file +11 0.60000001 +12 99