diff --git a/README.md b/README.md index ea76999..a718271 100644 --- a/README.md +++ b/README.md @@ -4,4 +4,4 @@ TODO - Error analysis - Better error handling (file read/write) - Better File I/O performance (file read in slow) -- Output to file +- integration tests diff --git a/src/cli.yml b/src/cli.yml index 7e241cb..2a12877 100644 --- a/src/cli.yml +++ b/src/cli.yml @@ -54,4 +54,10 @@ args: help: WHAM temperature in Kelvin takes_value: true required: true + - output: + short: o + long: output + help: Free energy output file (defaults to "wham.out") + takes_value: true + required: false diff --git a/src/io.rs b/src/io.rs index 83d756f..86ee952 100644 --- a/src/io.rs +++ b/src/io.rs @@ -8,6 +8,8 @@ use k_B; use std::process; use std::option::Option; use std::path::Path; +use std::error::Error; +use std::result::Result; // Returns the path to path2 relative to path1 // path1: "path/to/file.dat" @@ -133,6 +135,16 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option { } } +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() { + let x = ds.get_x_for_bin(bin); + writeln!(output, "{:8.6}\t{:8.6}\t{:8.6}", x, free[bin], prob[bin])?; + } + Ok(()) +} + #[cfg(test)] mod tests { use super::*; @@ -147,7 +159,8 @@ mod tests { tolerance: 0.0, max_iterations: 0, temperature: 300.0, - cyclic: false + cyclic: false, + output: "qwert".to_string(), } } diff --git a/src/lib.rs b/src/lib.rs index c33bde1..aed0691 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -27,13 +27,14 @@ pub struct Config { pub max_iterations: usize, pub temperature: f32, pub cyclic: bool, + pub output: String, } 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={}" , 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.max_iterations, self.temperature, self.cyclic) } } @@ -170,14 +171,15 @@ fn free_energy(ds: &Dataset, P: &mut Vec, A: &mut Vec) { for bin in 0..ds.num_bins { A[bin] -= bin_min; } + // Normalize P - // let mut P_sum = 0.0; - // for bin in 0..ds.num_bins { - // P_sum += P[bin]; - // } - // for bin in 0..ds.num_bins { - // P[bin] /= P_sum; - // } + let mut P_sum = 0.0; + for bin in 0..ds.num_bins { + P_sum += P[bin]; + } + for bin in 0..ds.num_bins { + P[bin] /= P_sum; + } } @@ -221,12 +223,15 @@ pub fn run(cfg: &Config) -> Result<(), Box>{ } // final free energy calculation and state dump + println!("Finished. Dumping final PMF"); free_energy(&histograms, &mut P, &mut A); dump_state(&histograms, &F, &F_prev, &P, &A); if iteration == cfg.max_iterations { println!("!!!!! WHAM not converged! (max iterations reached) !!!!!"); } + + io::write_results(&cfg.output, &histograms, &A, &P)?; Ok(()) } diff --git a/src/main.rs b/src/main.rs index cf6250e..7c486ce 100644 --- a/src/main.rs +++ b/src/main.rs @@ -21,9 +21,10 @@ fn cli() -> Result> { let tolerance: f32 = 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(); Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, - verbose, tolerance, max_iterations, temperature, cyclic}) + verbose, tolerance, max_iterations, temperature, cyclic, output}) } fn main() {