bauer/WHAM
1
0
Fork 0
mirror of https://github.com/dnlbauer/WHAM.git synced 2026-09-10 22:25:31 +00:00
Code Issues Packages Projects Releases Wiki Activity

partial io for multidimensional WHAM

Browse Source
This commit is contained in:
Daniel Bauer
2018-10-13 15:07:06 +02:00
parent 6450839c4a
commit 29d2dd2b03
4 changed files with 73 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/histogram.rs
View File

@@ -62,7 +62,7 @@ pub struct Dataset {
pub cyclic: bool, pub cyclic: bool,
// locations of biases // locations of biases
bias_x0: Vec<f64>, bias_pos: Vec<f64>,
// force constants of biases // force constants of biases
bias_fc: Vec<f64>, bias_fc: Vec<f64>,
@@ -73,7 +73,7 @@ pub struct Dataset {
impl Dataset { impl 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 { pub fn new(num_bins: usize, bin_width: f64, hist_min: f64, hist_max: f64, bias_pos: 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<f64>>> = RefCell::new(vec![None; num_bins*num_windows]); let bias: RefCell<Vec<Option<f64>>> = RefCell::new(vec![None; num_bins*num_windows]);
Dataset{ Dataset{
@@ -85,7 +85,7 @@ impl Dataset {
kT, kT,
histograms, histograms,
cyclic, cyclic,
bias_x0, bias_pos,
bias_fc, bias_fc,
bias, bias,
} }
@@ -102,7 +102,7 @@ impl Dataset {
Some(val) => val, Some(val) => val,
None => { None => {
let x = self.get_x_for_bin(bin); 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 { if self.cyclic {
let hist_len = self.hist_max-self.hist_min; let hist_len = self.hist_max-self.hist_min;
if dx > 0.5*hist_len { if dx > 0.5*hist_len {

71
src/io.rs
View File

@@ -29,7 +29,7 @@ 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<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 histograms: Vec<Histogram> = Vec::new(); let mut histograms: Vec<Histogram> = Vec::new();
let kT = cfg.temperature * k_B; let kT = cfg.temperature * k_B;
@@ -45,9 +45,11 @@ 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, 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) { match read_window_file(&path, cfg) {
Some(h) => { Some(h) => {
histograms.push(h); histograms.push(h);
@@ -55,25 +57,37 @@ pub fn read_data(cfg: &Config) -> Option<Dataset> {
}, },
None => { None => {
eprintln!("No data points inside histogram boundaries: {}", &path); 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(|| { // parse bias force constants and positions
eprintln!("Failed to read coordinate from: {}", &line); for _ in 0..cfg.dimens {
process::exit(1) match split.next()?.parse() {
})); Ok(x) => bias_pos.push(x),
bias_fc.push(k.unwrap_or_else(|| { _ => {
eprintln!("Failed to read bias value from: {}", &line); eprintln!("Failed to read bias coordinate.");
process::exit(1) 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 { if histograms.len() > 0 {
let bin_width = (cfg.hist_max - cfg.hist_min)/(cfg.num_bins as f64); let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
Some(Dataset::new(cfg.num_bins, bin_width, cfg.hist_min, cfg.hist_max, bias_x0, bias_fc, kT, histograms, cfg.cyclic)) (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 { } else {
None None
} }
@@ -87,8 +101,8 @@ 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[0]];
let bin_width = (cfg.hist_max - cfg.hist_min)/(cfg.num_bins as f64); let bin_width = (cfg.hist_max[0] - cfg.hist_min[0])/(cfg.num_bins[0] 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
@@ -100,8 +114,8 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option<Histogram> {
match x { match x {
Some(x) => { Some(x) => {
if x > cfg.hist_min && x < cfg.hist_max { if x > cfg.hist_min[0] && x < cfg.hist_max[0] {
let bin_ndx = ((x-cfg.hist_min) / bin_width) as usize; let bin_ndx = ((x-cfg.hist_min[0]) / bin_width) as usize;
global_hist[bin_ndx] += 1.0; global_hist[bin_ndx] += 1.0;
} }
} }
@@ -113,7 +127,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option<Histogram> {
} }
let mut max_bin: usize = 0; 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() { for bin in 0..global_hist.len() {
if global_hist[bin] != 0.0 && bin > max_bin { if global_hist[bin] != 0.0 && bin > max_bin {
max_bin = bin; max_bin = bin;
@@ -152,9 +166,10 @@ mod tests {
fn cfg() -> Config { fn cfg() -> Config {
Config { Config {
metadata_file: "tests/data/metadata.dat".to_string(), metadata_file: "tests/data/metadata.dat".to_string(),
hist_min: 0.0, hist_min: vec![0.0],
hist_max: 3.0, hist_max: vec![3.0],
num_bins: 30, num_bins: vec![30],
dimens: 1,
verbose: false, verbose: false,
tolerance: 0.0, tolerance: 0.0,
max_iterations: 0, max_iterations: 0,
@@ -187,13 +202,13 @@ mod tests {
let ds = ds.unwrap(); let ds = ds.unwrap();
println!("{:?}", ds); println!("{:?}", ds);
assert_eq!(2, ds.num_windows); 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 // fields are private
// assert_eq!(cfg.hist_min, ds.hist_min); // assert_eq!(cfg.hist_min[0], ds.hist_min[0]);
// assert_eq!(cfg.hist_max, ds.hist_max); // assert_eq!(cfg.hist_max[0], ds.hist_max[0]);
// let expected_bin_width = (cfg.hist_max - cfg.hist_min)/cfg.num_bins as f64; // 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!(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!(vec![100.0, 200.0], ds.bias_fc);
assert_eq!(cfg.temperature * k_B, ds.kT); assert_eq!(cfg.temperature * k_B, ds.kT);
assert_eq!(2, ds.histograms.len()) assert_eq!(2, ds.histograms.len())

13
src/lib.rs
View File

@@ -19,9 +19,10 @@ static k_B: f64 = 0.0083144621; // kJ/mol*K
#[derive(Debug)] #[derive(Debug)]
pub struct Config { pub struct Config {
pub metadata_file: String, pub metadata_file: String,
pub hist_min: f64, pub hist_min: Vec<f64>,
pub hist_max: f64, pub hist_max: Vec<f64>,
pub num_bins: usize, pub num_bins: Vec<usize>,
pub dimens: usize,
pub verbose: bool, pub verbose: bool,
pub tolerance: f64, pub tolerance: f64,
pub max_iterations: usize, pub max_iterations: usize,
@@ -32,7 +33,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={}\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.hist_max, self.num_bins, self.verbose, self.tolerance,
self.max_iterations, self.temperature, self.cyclic) self.max_iterations, self.temperature, self.cyclic)
} }
@@ -95,7 +96,7 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec<f64>,F: &mut Vec<f64>, P: &
// } // }
// let bias = calc_bias( // let bias = calc_bias(
// ds.bias_fc[window], // ds.bias_fc[window],
// ds.bias_x0[window], // ds.bias_pos[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 f64* bf // denom += ds.histograms[window].num_points as f64* bf
@@ -105,7 +106,7 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec<f64>,F: &mut Vec<f64>, P: &
// for window in 0..ds.num_windows { // for window in 0..ds.num_windows {
// let bias = calc_bias( // let bias = calc_bias(
// ds.bias_fc[window], // ds.bias_fc[window],
// ds.bias_x0[window], // ds.bias_pos[window],
// x); // x);
// let bf = (-bias/ds.kT).exp() * P[bin]; // let bf = (-bias/ds.kT).exp() * P[bin];
// F[window] += bf; // F[window] += bf;

22
src/main.rs
View File

@@ -13,17 +13,29 @@ 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: 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 verbose: bool = matches.is_present("verbose");
let temperature: f64 = matches.value_of("temperature").unwrap().parse()?; let temperature: f64 = matches.value_of("temperature").unwrap().parse()?;
let tolerance: f64 = 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 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");
Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, let hist_min: Vec<f64> = matches.value_of("min_hist").unwrap()
.split(':').map(|x| { x.parse().unwrap() }).collect();
let hist_max: Vec<f64> = matches.value_of("max_hist").unwrap()
.split(':').map(|x| { x.parse().unwrap() }).collect();
let num_bins: Vec<usize> = 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}) verbose, tolerance, max_iterations, temperature, cyclic, output})
} }