diff --git a/src/histogram.rs b/src/histogram.rs index 0c0ef66..cd4f02f 100644 --- a/src/histogram.rs +++ b/src/histogram.rs @@ -62,7 +62,7 @@ pub struct Dataset { pub cyclic: bool, // locations of biases - bias_x0: Vec, + bias_pos: Vec, // force constants of biases bias_fc: Vec, @@ -73,7 +73,7 @@ pub struct Dataset { impl 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 { + pub fn new(num_bins: usize, bin_width: f64, hist_min: f64, hist_max: f64, bias_pos: 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]); Dataset{ @@ -85,7 +85,7 @@ impl Dataset { kT, histograms, cyclic, - bias_x0, + bias_pos, bias_fc, bias, } @@ -102,7 +102,7 @@ impl Dataset { Some(val) => val, None => { let x = self.get_x_for_bin(bin); - let mut dx = (x-self.bias_x0[window]).abs(); + let mut dx = (x-self.bias_pos[window]).abs(); if self.cyclic { let hist_len = self.hist_max-self.hist_min; if dx > 0.5*hist_len { diff --git a/src/io.rs b/src/io.rs index 432a397..bf01a57 100644 --- a/src/io.rs +++ b/src/io.rs @@ -29,7 +29,7 @@ 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_pos: Vec = Vec::new(); let mut bias_fc: Vec = Vec::new(); let mut histograms: Vec = Vec::new(); let kT = cfg.temperature * k_B; @@ -45,9 +45,11 @@ pub fn read_data(cfg: &Config) -> Option { if line.starts_with("#") || line.len() == 0 { continue; } - let (path, x0, k) = scan_fmt!(&line, "{} {} {}", String, f64, f64); - - let path = get_relative_path(&cfg.metadata_file, &path.unwrap()); + + let mut split = line.split_whitespace(); + + // parse histogram data + let path = get_relative_path(&cfg.metadata_file, split.next()?); match read_window_file(&path, cfg) { Some(h) => { histograms.push(h); @@ -55,25 +57,37 @@ pub fn read_data(cfg: &Config) -> Option { }, None => { eprintln!("No data points inside histogram boundaries: {}", &path); - continue; // goto next iteration and skip adding bias values + process::exit(1) } } - bias_x0.push(x0.unwrap_or_else(|| { - eprintln!("Failed to read coordinate from: {}", &line); - process::exit(1) - })); - bias_fc.push(k.unwrap_or_else(|| { - eprintln!("Failed to read bias value from: {}", &line); - process::exit(1) - })); + // parse bias force constants and positions + for _ in 0..cfg.dimens { + match split.next()?.parse() { + Ok(x) => bias_pos.push(x), + _ => { + eprintln!("Failed to read bias coordinate."); + process::exit(1); + } + } + } + for _ in 0..cfg.dimens { + match split.next()?.parse() { + Ok(x) => bias_fc.push(x), + _ => { + eprintln!("Failed to read bias force constant."); + process::exit(1); + } + } - + } } if histograms.len() > 0 { - 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)) + let bin_width: Vec = (0..cfg.dimens).map(|idx| { + (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64) + }).collect(); + Some(Dataset::new(cfg.num_bins[0], bin_width[0], cfg.hist_min[0], cfg.hist_max[0], bias_pos, bias_fc, kT, histograms, cfg.cyclic)) } else { None } @@ -87,8 +101,8 @@ 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 f64); + let mut global_hist = vec![0.0; cfg.num_bins[0]]; + let bin_width = (cfg.hist_max[0] - cfg.hist_min[0])/(cfg.num_bins[0] as f64); for l in buf.lines() { let line = l.unwrap(); // skip comments and empty lines @@ -100,8 +114,8 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option { match x { Some(x) => { - if x > cfg.hist_min && x < cfg.hist_max { - let bin_ndx = ((x-cfg.hist_min) / bin_width) as usize; + if x > cfg.hist_min[0] && x < cfg.hist_max[0] { + let bin_ndx = ((x-cfg.hist_min[0]) / bin_width) as usize; global_hist[bin_ndx] += 1.0; } } @@ -113,7 +127,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option { } let mut max_bin: usize = 0; - let mut min_bin: usize =(cfg.num_bins-1) as usize; + let mut min_bin: usize =(cfg.num_bins[0]-1) as usize; for bin in 0..global_hist.len() { if global_hist[bin] != 0.0 && bin > max_bin { max_bin = bin; @@ -150,11 +164,12 @@ mod tests { use super::*; fn cfg() -> Config { - Config{ + Config { metadata_file: "tests/data/metadata.dat".to_string(), - hist_min: 0.0, - hist_max: 3.0, - num_bins: 30, + hist_min: vec![0.0], + hist_max: vec![3.0], + num_bins: vec![30], + dimens: 1, verbose: false, tolerance: 0.0, max_iterations: 0, @@ -187,13 +202,13 @@ mod tests { let ds = ds.unwrap(); println!("{:?}", ds); assert_eq!(2, ds.num_windows); - assert_eq!(cfg.num_bins, ds.num_bins); + assert_eq!(cfg.num_bins[0], ds.num_bins[0]); // fields are private - // 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 f64; + // assert_eq!(cfg.hist_min[0], ds.hist_min[0]); + // assert_eq!(cfg.hist_max[0], ds.hist_max[0]); + // let expected_bin_width = (cfg.hist_max[0] - cfg.hist_min[0])/cfg.num_bins[0] as f64; // assert_eq!(expected_bin_width, ds.bin_width); - // assert_eq!(vec![0.0, 1.0], ds.bias_x0); + // assert_eq!(vec![0.0, 1.0], ds.bias_pos); // assert_eq!(vec![100.0, 200.0], ds.bias_fc); assert_eq!(cfg.temperature * k_B, ds.kT); assert_eq!(2, ds.histograms.len()) diff --git a/src/lib.rs b/src/lib.rs index ebc7b61..d52d792 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -19,9 +19,10 @@ static k_B: f64 = 0.0083144621; // kJ/mol*K #[derive(Debug)] pub struct Config { pub metadata_file: String, - pub hist_min: f64, - pub hist_max: f64, - pub num_bins: usize, + pub hist_min: Vec, + pub hist_max: Vec, + pub num_bins: Vec, + pub dimens: usize, pub verbose: bool, pub tolerance: f64, pub max_iterations: usize, @@ -32,7 +33,7 @@ pub struct Config { impl fmt::Display for Config { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!(f, "Metadata={}, hist_min={}, hist_max={}, bins={}\nverbose={}, tolerance={}, iterations={}, temperature={}, cyclic={}" , self.metadata_file, self.hist_min, + write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?}\nverbose={}, tolerance={}, iterations={}, temperature={}, cyclic={:?}" , self.metadata_file, self.hist_min, self.hist_max, self.num_bins, self.verbose, self.tolerance, self.max_iterations, self.temperature, self.cyclic) } @@ -95,7 +96,7 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec,F: &mut Vec, P: & // } // let bias = calc_bias( // ds.bias_fc[window], - // ds.bias_x0[window], + // ds.bias_pos[window], // x); // let bf = ((F_prev[window]-bias) / ds.kT).exp(); // denom += ds.histograms[window].num_points as f64* bf @@ -105,7 +106,7 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec,F: &mut Vec, P: & // for window in 0..ds.num_windows { // let bias = calc_bias( // ds.bias_fc[window], - // ds.bias_x0[window], + // ds.bias_pos[window], // x); // let bf = (-bias/ds.kT).exp() * P[bin]; // F[window] += bf; diff --git a/src/main.rs b/src/main.rs index c4b1844..a5b95f7 100644 --- a/src/main.rs +++ b/src/main.rs @@ -13,17 +13,29 @@ 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: 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: 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(); + let cyclic: bool = matches.is_present("cyclic"); - Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, + let hist_min: Vec = matches.value_of("min_hist").unwrap() + .split(':').map(|x| { x.parse().unwrap() }).collect(); + let hist_max: Vec = matches.value_of("max_hist").unwrap() + .split(':').map(|x| { x.parse().unwrap() }).collect(); + let num_bins: Vec = matches.value_of("bins").unwrap() + .split(':').map(|x| { x.parse().unwrap() }).collect(); + + if num_bins.len() != hist_max.len() || num_bins.len() != hist_max.len() { + eprintln!("Input dimensions do not match (min: {}, max: {}, bins: {})", + hist_min.len(), hist_max.len(), num_bins.len()); + process::exit(1); + } + + let dimens = num_bins.len(); + + Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens, verbose, tolerance, max_iterations, temperature, cyclic, output}) }