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

use 64 bit floats

Browse Source
This commit is contained in:
Daniel Bauer
2018-10-08 10:13:19 +02:00
parent 23085ab6f1
commit 88d9146911
5 changed files with 67 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files

34
src/histogram.rs
View File

@@ -14,18 +14,18 @@ pub struct Histogram {
pub num_points: u32, pub num_points: u32,
// histogram bins // histogram bins
pub bins: Vec<f32> pub bins: Vec<f64>
} }
impl Histogram { impl Histogram {
pub fn new(first: usize, last: usize, num_points: u32, bins: Vec<f32>) -> Histogram { pub fn new(first: usize, last: usize, num_points: u32, bins: Vec<f64>) -> Histogram {
assert_eq!(last-first+1, bins.len(), "histogram length does not match first/last."); assert_eq!(last-first+1, bins.len(), "histogram length does not match first/last.");
Histogram {first, last, num_points, bins} Histogram {first, last, num_points, bins}
} }
// Returns the value of a bin if the bin is present in this // Returns the value of a bin if the bin is present in this
// histogram // histogram
pub fn get_bin_count(&self, bin: usize) -> Option<f32> { pub fn get_bin_count(&self, bin: usize) -> Option<f64> {
if bin < self.first || bin > self.last { if bin < self.first || bin > self.last {
None None
} else { } else {
@@ -44,16 +44,16 @@ pub struct Dataset {
pub num_bins: usize, pub num_bins: usize,
// min value of the histogram // min value of the histogram
pub hist_min: f32, pub hist_min: f64,
// max value of the histogram // max value of the histogram
pub hist_max: f32, pub hist_max: f64,
// width of a bin in unit of x // width of a bin in unit of x
pub bin_width: f32, pub bin_width: f64,
// value of kT // value of kT
pub kT: f32, pub kT: f64,
// histogram for each window // histogram for each window
pub histograms: Vec<Histogram>, pub histograms: Vec<Histogram>,
@@ -62,20 +62,20 @@ pub struct Dataset {
pub cyclic: bool, pub cyclic: bool,
// locations of biases // locations of biases
bias_x0: Vec<f32>, bias_x0: Vec<f64>,
// force constants of biases // force constants of biases
bias_fc: Vec<f32>, bias_fc: Vec<f64>,
// bias value cache // bias value cache
bias: RefCell<Vec<Option<f32>>> bias: RefCell<Vec<Option<f64>>>
} }
impl Dataset { impl Dataset {
pub fn new(num_bins: usize, bin_width: f32, hist_min: f32, hist_max: f32, bias_x0: Vec<f32>, bias_fc: Vec<f32>, kT: f32, histograms: Vec<Histogram>, cyclic: bool) -> Dataset { pub fn new(num_bins: usize, bin_width: f64, hist_min: f64, hist_max: f64, bias_x0: Vec<f64>, bias_fc: Vec<f64>, kT: f64, histograms: Vec<Histogram>, cyclic: bool) -> Dataset {
let num_windows = histograms.len(); let num_windows = histograms.len();
let bias: RefCell<Vec<Option<f32>>> = RefCell::new(vec![None; num_bins*num_windows]); let bias: RefCell<Vec<Option<f64>>> = RefCell::new(vec![None; num_bins*num_windows]);
Dataset{ Dataset{
num_windows, num_windows,
num_bins, num_bins,
@@ -95,7 +95,7 @@ impl Dataset {
// 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 // 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 ndx = bin + (self.num_bins*window);
let mut cache = self.bias.borrow_mut(); let mut cache = self.bias.borrow_mut();
match cache[ndx] { match cache[ndx] {
@@ -117,8 +117,8 @@ impl Dataset {
} }
// get center x value for a bin // get center x value for a bin
pub fn get_x_for_bin(&self, bin: usize) -> f32 { pub fn get_x_for_bin(&self, bin: usize) -> f64 {
self.hist_min + self.bin_width * ((bin as f32) + 0.5) self.hist_min + self.bin_width * ((bin as f64) + 0.5)
} }
} }
@@ -212,8 +212,8 @@ mod tests {
#[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<f32> = 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 f32 + 0.5).collect(); .map(|x| *x as f64 + 0.5).collect();
for i in 0..9 { for i in 0..9 {
assert_eq!(expected[i], ds.get_x_for_bin(i)); assert_eq!(expected[i], ds.get_x_for_bin(i));
} }

18
src/io.rs
View File

@@ -29,8 +29,8 @@ pub fn vprintln(s: String, verbose: bool) {
// Read input data into a histogram set by iterating over input files // Read input data into a histogram set by iterating over input files
// given in the metadata file // given in the metadata file
pub fn read_data(cfg: &Config) -> Option<Dataset> { pub fn read_data(cfg: &Config) -> Option<Dataset> {
let mut bias_x0: Vec<f32> = Vec::new(); let mut bias_x0: Vec<f64> = Vec::new();
let mut bias_fc: Vec<f32> = Vec::new(); let mut bias_fc: Vec<f64> = Vec::new();
let mut histograms: Vec<Histogram> = Vec::new(); let mut histograms: Vec<Histogram> = Vec::new();
let kT = cfg.temperature * k_B; let kT = cfg.temperature * k_B;
let f = File::open(&cfg.metadata_file).unwrap_or_else(|x| { let f = File::open(&cfg.metadata_file).unwrap_or_else(|x| {
@@ -45,7 +45,7 @@ pub fn read_data(cfg: &Config) -> Option<Dataset> {
if line.starts_with("#") || line.len() == 0 { if line.starts_with("#") || line.len() == 0 {
continue; 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()); let path = get_relative_path(&cfg.metadata_file, &path.unwrap());
match read_window_file(&path, cfg) { match read_window_file(&path, cfg) {
@@ -72,7 +72,7 @@ pub fn read_data(cfg: &Config) -> Option<Dataset> {
} }
if histograms.len() > 0 { 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)) Some(Dataset::new(cfg.num_bins, bin_width, cfg.hist_min, cfg.hist_max, bias_x0, bias_fc, kT, histograms, cfg.cyclic))
} else { } else {
None None
@@ -88,7 +88,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option<Histogram> {
let buf = BufReader::new(&f); let buf = BufReader::new(&f);
let mut global_hist = vec![0.0; cfg.num_bins]; 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() { for l in buf.lines() {
let line = l.unwrap(); let line = l.unwrap();
// skip comments and empty lines // skip comments and empty lines
@@ -96,7 +96,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option<Histogram> {
continue; continue;
} }
let (_, x) = scan_fmt!(&line, "{} {}", f32, f32); let (_, x) = scan_fmt!(&line, "{} {}", f64, f64);
match x { match x {
Some(x) => { Some(x) => {
@@ -130,12 +130,12 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option<Histogram> {
global_hist.truncate(max_bin+1); global_hist.truncate(max_bin+1);
global_hist.drain(..min_bin); 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)) Some(Histogram::new(min_bin, max_bin, num_points as u32, global_hist))
} }
} }
pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec<f32>, prob: &Vec<f32>) -> Result<(), Box<Error>> { pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec<f64>, prob: &Vec<f64>) -> Result<(), Box<Error>> {
let mut output = File::create(out_file)?; let mut output = File::create(out_file)?;
writeln!(output, "#{:8}\t{:8}\t{:8}", "x", "Free Energy", "Probability"); writeln!(output, "#{:8}\t{:8}\t{:8}", "x", "Free Energy", "Probability");
for bin in 0..free.len() { for bin in 0..free.len() {
@@ -191,7 +191,7 @@ mod tests {
assert_eq!(cfg.num_bins, ds.num_bins); assert_eq!(cfg.num_bins, ds.num_bins);
assert_eq!(cfg.hist_min, ds.hist_min); assert_eq!(cfg.hist_min, ds.hist_min);
assert_eq!(cfg.hist_max, ds.hist_max); 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); assert_eq!(expected_bin_width, ds.bin_width);
// bias fields are private // bias fields are private
// assert_eq!(vec![0.0, 1.0], ds.bias_x0); // assert_eq!(vec![0.0, 1.0], ds.bias_x0);

60
src/lib.rs
View File

@@ -9,23 +9,23 @@ pub mod histogram;
use std::error::Error; use std::error::Error;
use std::result::Result; use std::result::Result;
use histogram::{Dataset,Histogram}; use histogram::{Dataset,Histogram};
use std::f32; use std::f64;
use std::fmt; use std::fmt;
#[allow(non_upper_case_globals)] #[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 // Application config
#[derive(Debug)] #[derive(Debug)]
pub struct Config { pub struct Config {
pub metadata_file: String, pub metadata_file: String,
pub hist_min: f32, pub hist_min: f64,
pub hist_max: f32, pub hist_max: f64,
pub num_bins: usize, pub num_bins: usize,
pub verbose: bool, pub verbose: bool,
pub tolerance: f32, pub tolerance: f64,
pub max_iterations: usize, pub max_iterations: usize,
pub temperature: f32, pub temperature: f64,
pub cyclic: bool, pub cyclic: bool,
pub output: String, pub output: String,
} }
@@ -41,7 +41,7 @@ impl fmt::Display for Config {
// Checks for convergence between two WHAM iterations. WHAM is considered as // Checks for convergence between two WHAM iterations. WHAM is considered as
// converged if the absolute difference for the calculated bias offset is // converged if the absolute difference for the calculated bias offset is
// smaller then a tolerance value for every simulation window. // smaller then a tolerance value for every simulation window.
fn is_converged(old_F: &Vec<f32>, new_F: &Vec<f32>, tolerance: f32) -> bool { fn is_converged(old_F: &Vec<f64>, new_F: &Vec<f64>, tolerance: f64) -> bool {
!new_F.iter().zip(old_F.iter()) !new_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 })
@@ -49,9 +49,9 @@ fn is_converged(old_F: &Vec<f32>, new_F: &Vec<f32>, tolerance: f32) -> bool {
// estimate the probability of a bin of the histogram set based on F values // estimate the probability of a bin of the histogram set based on F values
// This evaluates the first WHAM equation for each bin // This evaluates the first WHAM equation for each bin
fn calc_bin_probability(bin: usize, ds: &Dataset, F: &Vec<f32>) -> f32 { fn calc_bin_probability(bin: usize, ds: &Dataset, F: &Vec<f64>) -> f64 {
let mut denom_sum = 0.0; let mut denom_sum: f64 = 0.0;
let mut bin_count = 0.0; let mut bin_count: f64 = 0.0;
for window in 0..ds.num_windows { for window in 0..ds.num_windows {
let h: &Histogram = &ds.histograms[window]; let h: &Histogram = &ds.histograms[window];
if let Some(count) = h.get_bin_count(bin) { if let Some(count) = h.get_bin_count(bin) {
@@ -59,16 +59,16 @@ fn calc_bin_probability(bin: usize, ds: &Dataset, F: &Vec<f32>) -> f32 {
} }
let bias = ds.calc_bias(bin, window); let bias = ds.calc_bias(bin, window);
let bias_offset = ((F[window] - bias) / ds.kT).exp(); 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 bin_count / denom_sum
} }
// estimate the bias offset F of the histogram based on given probabilities // estimate the bias offset F of the histogram based on given probabilities
// This evaluates the second WHAM equation for each window // This evaluates the second WHAM equation for each window
fn calc_window_F(window: usize, ds: &Dataset, P: &Vec<f32>) -> f32 { fn calc_window_F(window: usize, ds: &Dataset, P: &Vec<f64>) -> f64 {
let bf_sum: f32 = (0..ds.num_bins).zip(P.iter()) // zip bins and P let bf_sum: f64 = (0..ds.num_bins).zip(P.iter()) // zip bins and P
.map(|x: (usize, &f32)| { .map(|x: (usize, &f64)| {
x.1 * (-ds.calc_bias(x.0, window)/ds.kT).exp() x.1 * (-ds.calc_bias(x.0, window)/ds.kT).exp()
}).sum(); }).sum();
-ds.kT * bf_sum.ln() -ds.kT * bf_sum.ln()
@@ -77,7 +77,7 @@ fn calc_window_F(window: usize, ds: &Dataset, P: &Vec<f32>) -> f32 {
// One full WHAM iteration includes calculation of new probabilities P and // One full WHAM iteration includes calculation of new probabilities P and
// new bias offsets F based on previous bias offsets F_prev. This updates // new bias offsets F based on previous bias offsets F_prev. This updates
// the values in vectors F and P // the values in vectors F and P
fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec<f32>,F: &mut Vec<f32>, P: &mut Vec<f32>) { fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec<f64>,F: &mut Vec<f64>, P: &mut Vec<f64>) {
// reset bias offsets // reset bias offsets
for window in 0..ds.num_windows { for window in 0..ds.num_windows {
F[window] = 0.0; F[window] = 0.0;
@@ -98,7 +98,7 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec<f32>,F: &mut Vec<f32>, P: &
// ds.bias_x0[window], // ds.bias_x0[window],
// x); // x);
// let bf = ((F_prev[window]-bias) / ds.kT).exp(); // 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; // P[bin] = num / denom;
@@ -141,18 +141,18 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec<f32>,F: &mut Vec<f32>, P: &
} }
// get average difference between two bias offset sets // get average difference between two bias offset sets
fn diff_avg(F: &Vec<f32>, F_prev: &Vec<f32>) -> f32 { fn diff_avg(F: &Vec<f64>, F_prev: &Vec<f64>) -> f64 {
let mut F_sum = 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 f32 F_sum / F.len() as f64
} }
// calculate the normalized free energy from normalized probability values // calculate the normalized free energy from normalized probability values
fn free_energy(ds: &Dataset, P: &mut Vec<f32>, A: &mut Vec<f32>) { fn free_energy(ds: &Dataset, P: &mut Vec<f64>, A: &mut Vec<f64>) {
let mut bin_min = f32::MAX; let mut bin_min = f64::MAX;
// Free energy calculation // Free energy calculation
for bin in 0..ds.num_bins { for bin in 0..ds.num_bins {
@@ -190,10 +190,10 @@ pub fn run(cfg: &Config) -> Result<(), Box<Error>>{
println!("{}",&histograms); println!("{}",&histograms);
// allocate only once for better performance // allocate only once for better performance
let mut F_prev = vec![f32::INFINITY; histograms.num_windows]; let mut F_prev: Vec<f64> = vec![f64::INFINITY; histograms.num_windows];
let mut F = vec![0.0; histograms.num_windows]; let mut F: Vec<f64> = vec![0.0; histograms.num_windows];
let mut P = vec![f32::NAN; histograms.num_bins]; let mut P: Vec<f64> = vec![f64::NAN; histograms.num_bins];
let mut A = vec![f32::NAN; histograms.num_bins]; let mut A: Vec<f64> = vec![f64::NAN; histograms.num_bins];
// perform WHAM until convergence // perform WHAM until convergence
let mut iteration = 0; let mut iteration = 0;
@@ -240,7 +240,7 @@ pub fn run(cfg: &Config) -> Result<(), Box<Error>>{
Ok(()) Ok(())
} }
fn dump_state(ds: &Dataset, F: &Vec<f32>, F_prev: &Vec<f32>, P: &Vec<f32>, A: &Vec<f32>) { fn dump_state(ds: &Dataset, F: &Vec<f64>, F_prev: &Vec<f64>, P: &Vec<f64>, A: &Vec<f64>) {
println!("# PMF"); println!("# PMF");
println!("#x\t\tFree Energy\t\tP(x)"); println!("#x\t\tFree Energy\t\tP(x)");
for bin in 0..ds.num_bins { for bin in 0..ds.num_bins {
@@ -257,7 +257,7 @@ fn dump_state(ds: &Dataset, F: &Vec<f32>, F_prev: &Vec<f32>, P: &Vec<f32>, A: &V
#[cfg(test)] #[cfg(test)]
mod tests { mod tests {
use super::histogram::{Dataset,Histogram}; use super::histogram::{Dataset,Histogram};
use std::f32; use std::f64;
#[test] #[test]
fn is_converged() { 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) 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(); let d = (a-b).abs();
assert!(d <= tolerance, "Values are not close: {}, {}, d={}", &a, &b, &d); assert!(d <= tolerance, "Values are not close: {}, {}, d={}", &a, &b, &d);
} }
@@ -317,9 +317,9 @@ mod tests {
let ds = create_test_ds(); let ds = create_test_ds();
let prev_F = vec![0.0; ds.num_windows]; let prev_F = vec![0.0; ds.num_windows];
let mut 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); 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); let expected_P = vec!(0.959, 0.331, 0.656, 46.750);
for bin in 0..ds.num_bins { for bin in 0..ds.num_bins {
assert_near(expected_P[bin], P[bin], 0.01) assert_near(expected_P[bin], P[bin], 0.01)

8
src/main.rs
View File

@@ -13,12 +13,12 @@ fn cli() -> Result<Config, Box<Error>> {
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 hist_min: f32 = matches.value_of("min_hist").unwrap().parse()?; let hist_min: f64 = matches.value_of("min_hist").unwrap().parse()?;
let hist_max: f32 = matches.value_of("max_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 num_bins: usize = matches.value_of("bins").unwrap().parse()?;
let verbose: bool = matches.is_present("verbose"); let verbose: bool = matches.is_present("verbose");
let temperature: f32 = matches.value_of("temperature").unwrap().parse()?; let temperature: f64 = matches.value_of("temperature").unwrap().parse()?;
let tolerance: f32 = matches.value_of("tolerance").unwrap_or("0.000001").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 max_iterations: usize = matches.value_of("iterations").unwrap_or("100000").parse()?;
let cyclic: bool = matches.is_present("cyclic"); let cyclic: bool = matches.is_present("cyclic");
let output = matches.value_of("output").unwrap_or("wham.out").to_string(); let output = matches.value_of("output").unwrap_or("wham.out").to_string();

12
tests/data/window_0.0.dat
View File

@@ -1,12 +1,12 @@
1 0.1 1 0.100000001
2 0.15 2 0.15
3 0.2 3 0.200000001
4 0.25 4 0.25
5 0.3 5 0.300000001
6 0.35 6 0.35
7 0.4 7 0.400000001
8 0.45 8 0.45
9 0.5 9 0.500000001
1 0.55 1 0.55
11 0.6 11 0.60000001
12 99 12 99