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

Merge branch 'precalc'

Browse Source
This commit is contained in:
Daniel Bauer
2018-10-18 01:03:39 +02:00
parent 7fd78fc2a5 c835520bf4
commit a13ebded6d
5 changed files with 10247 additions and 10273 deletions
Download Patch File Download Diff File Expand all files Collapse all files

166
example/1d/wham.out
View File

@@ -1,101 +1,101 @@
#x Free Energy Probability #x Free Energy Probability
-3.108600 7.118164 0.003561 -3.108600 7.118166 0.003561
-3.045800 5.331288 0.007290 -3.045800 5.331290 0.007290
-2.983000 3.879239 0.013048 -2.983000 3.879241 0.013048
-2.920200 2.968658 0.018796 -2.920200 2.968659 0.018796
-2.857400 1.942493 0.028362 -2.857400 1.942494 0.028362
-2.794600 1.418407 0.034994 -2.794600 1.418408 0.034994
-2.731800 1.169397 0.038668 -2.731800 1.169398 0.038668
-2.669000 0.843023 0.044073 -2.669000 0.843024 0.044073
-2.606200 0.636020 0.047887 -2.606200 0.636021 0.047887
-2.543400 0.720123 0.046299 -2.543400 0.720123 0.046299
-2.480600 1.066434 0.040297 -2.480600 1.066434 0.040297
-2.417800 1.509855 0.033734 -2.417800 1.509855 0.033734
-2.355000 2.015549 0.027544 -2.355000 2.015549 0.027544
-2.292200 2.226943 0.025306 -2.292200 2.226942 0.025306
-2.229400 2.573176 0.022026 -2.229400 2.573176 0.022026
-2.166600 2.591946 0.021861 -2.166600 2.591946 0.021861
-2.103800 2.476836 0.022893 -2.103800 2.476835 0.022893
-2.041000 2.513953 0.022555 -2.041000 2.513953 0.022555
-1.978200 2.470033 0.022956 -1.978200 2.470032 0.022956
-1.915400 2.227585 0.025299 -1.915400 2.227585 0.025299
-1.852600 2.105310 0.026570 -1.852600 2.105310 0.026570
-1.789800 1.778780 0.030286 -1.789800 1.778779 0.030286
-1.727000 1.447127 0.034593 -1.727000 1.447127 0.034593
-1.664200 0.906384 0.042968 -1.664200 0.906384 0.042968
-1.601400 0.304263 0.054699 -1.601400 0.304263 0.054699
-1.538600 0.247397 0.055960 -1.538600 0.247397 0.055960
-1.475800 0.000000 0.061795 -1.475800 0.000000 0.061795
-1.413000 0.536141 0.049843 -1.413000 0.536142 0.049843
-1.350200 1.428046 0.034859 -1.350200 1.428046 0.034859
-1.287400 2.398149 0.023627 -1.287400 2.398149 0.023627
-1.224600 3.784647 0.013552 -1.224600 3.784647 0.013552
-1.161800 5.696582 0.006297 -1.161800 5.696583 0.006297
-1.099000 7.656157 0.002870 -1.099000 7.656158 0.002870
-1.036200 9.954999 0.001142 -1.036200 9.955000 0.001142
-0.973400 12.357406 0.000436 -0.973400 12.357407 0.000436
-0.910600 14.954253 0.000154 -0.910600 14.954254 0.000154
-0.847800 17.745325 0.000050 -0.847800 17.745327 0.000050
-0.785000 20.559459 0.000016 -0.785000 20.559461 0.000016
-0.722200 22.783557 0.000007 -0.722200 22.783560 0.000007
-0.659400 25.175862 0.000003 -0.659400 25.175865 0.000003
-0.596600 26.419977 0.000002 -0.596600 26.419980 0.000002
-0.533800 27.885996 0.000001 -0.533800 27.885999 0.000001
-0.471000 29.029330 0.000001 -0.471000 29.029334 0.000001
-0.408200 30.366552 0.000000 -0.408200 30.366556 0.000000
-0.345400 31.615248 0.000000 -0.345400 31.615252 0.000000
-0.282600 32.781711 0.000000 -0.282600 32.781716 0.000000
-0.219800 33.731373 0.000000 -0.219800 33.731378 0.000000
-0.157000 34.476008 0.000000 -0.157000 34.476014 0.000000
-0.094200 35.387294 0.000000 -0.094200 35.387300 0.000000
-0.031400 35.579750 0.000000 -0.031400 35.579756 0.000000
0.031400 35.520188 0.000000 0.031400 35.520194 0.000000
0.094200 35.344757 0.000000 0.094200 35.344763 0.000000
0.157000 34.886453 0.000000 0.157000 34.886460 0.000000
0.219800 33.633992 0.000000 0.219800 33.633998 0.000000
0.282600 32.686504 0.000000 0.282600 32.686511 0.000000
0.345400 31.211770 0.000000 0.345400 31.211777 0.000000
0.408200 29.680443 0.000000 0.408200 29.680451 0.000000
0.471000 28.032233 0.000001 0.471000 28.032242 0.000001
0.533800 26.439591 0.000002 0.533800 26.439599 0.000002
0.596600 24.432458 0.000003 0.596600 24.432467 0.000003
0.659400 22.309884 0.000008 0.659400 22.309893 0.000008
0.722200 20.186031 0.000019 0.722200 20.186040 0.000019
0.785000 18.295037 0.000040 0.785000 18.295046 0.000040
0.847800 16.221223 0.000093 0.847800 16.221232 0.000093
0.910600 14.251901 0.000204 0.910600 14.251911 0.000204
0.973400 12.562343 0.000402 0.973400 12.562353 0.000402
1.036200 11.184019 0.000698 1.036200 11.184029 0.000698
1.099000 10.040667 0.001103 1.099000 10.040677 0.001103
1.161800 9.369054 0.001444 1.161800 9.369064 0.001444
1.224600 9.095382 0.001612 1.224600 9.095392 0.001612
1.287400 9.256025 0.001511 1.287400 9.256035 0.001511
1.350200 9.814058 0.001208 1.350200 9.814068 0.001208
1.413000 10.945197 0.000768 1.413000 10.945207 0.000768
1.475800 12.518888 0.000409 1.475800 12.518898 0.000409
1.538600 14.409735 0.000191 1.538600 14.409744 0.000191
1.601400 16.451956 0.000084 1.601400 16.451965 0.000084
1.664200 18.589157 0.000036 1.664200 18.589166 0.000036
1.727000 20.669496 0.000016 1.727000 20.669504 0.000016
1.789800 22.804783 0.000007 1.789800 22.804791 0.000007
1.852600 24.531682 0.000003 1.852600 24.531690 0.000003
1.915400 26.182604 0.000002 1.915400 26.182612 0.000002
1.978200 27.279432 0.000001 1.978200 27.279439 0.000001
2.041000 28.571373 0.000001 2.041000 28.571380 0.000001
2.103800 29.328963 0.000000 2.103800 29.328969 0.000000
2.166600 29.897732 0.000000 2.166600 29.897739 0.000000
2.229400 30.320901 0.000000 2.229400 30.320908 0.000000
2.292200 30.065385 0.000000 2.292200 30.065391 0.000000
2.355000 29.809130 0.000000 2.355000 29.809137 0.000000
2.417800 29.340595 0.000000 2.417800 29.340601 0.000000
2.480600 28.473341 0.000001 2.480600 28.473346 0.000001
2.543400 27.679348 0.000001 2.543400 27.679353 0.000001
2.606200 26.453155 0.000002 2.606200 26.453160 0.000002
2.669000 24.424644 0.000003 2.669000 24.424649 0.000003
2.731800 22.285929 0.000008 2.731800 22.285933 0.000008
2.794600 19.958403 0.000021 2.794600 19.958407 0.000021
2.857400 17.606287 0.000053 2.857400 17.606291 0.000053
2.920200 15.420135 0.000128 2.920200 15.420138 0.000128
2.983000 13.069596 0.000328 2.983000 13.069599 0.000328
3.045800 11.036083 0.000740 3.045800 11.036085 0.000740
3.108600 9.015655 0.001664 3.108600 9.015657 0.001664

20000
example/2d/wham.out
View File

File diff suppressed because it is too large Load Diff

29
src/histogram.rs
View File

@@ -126,6 +126,7 @@ impl Dataset {
} }
bias_sum += 0.5 * bias_fc[i] * dist * dist bias_sum += 0.5 * bias_fc[i] * dist * dist
} }
let bias_sum = (-bias_sum/self.kT).exp();
cache[ndx] = Some(bias_sum); cache[ndx] = Some(bias_sum);
bias_sum bias_sum
} }
@@ -149,6 +150,12 @@ mod tests {
use super::*; use super::*;
use super::super::k_B; use super::super::k_B;
macro_rules! assert_delta {
($x:expr, $y:expr, $d:expr) => {
assert!(($x-$y).abs() < $d, "{} != {}", $x, $y)
}
}
fn build_hist() -> Histogram { fn build_hist() -> Histogram {
Histogram::new( Histogram::new(
22, // num_points 22, // num_points
@@ -158,7 +165,7 @@ mod tests {
fn build_hist_set() -> Dataset { fn build_hist_set() -> Dataset {
let h = build_hist(); let h = build_hist();
Dataset::new( Dataset::new(
9, // num bins 9, // num bins
vec![1], vec![1],
vec![1.0], // bin width vec![1.0], // bin width
@@ -186,38 +193,40 @@ mod tests {
} }
} }
#[test] #[test]
fn calc_bias() { fn calc_bias() {
let ds = build_hist_set(); let ds = build_hist_set(); // k = 10
// 7th element -> x=7.5, x0=7.5 // 7th element -> x=7.5, x0=7.5
assert_eq!(0.0, ds.calc_bias(7, 0)); assert_delta!(1.0, ds.calc_bias(7, 0), 0.00000001);
// 8th element -> x=8.5, x0=7.5 // 8th element -> x=8.5, x0=7.5
assert_eq!(5.0, ds.calc_bias(8, 0)); assert_delta!(0.13472233779, ds.calc_bias(8,0), 0.00000001);
// 1st element -> x=0.5, x0=7.5. non-cyclic! // 1st element -> x=0.5, x0=7.5. non-cyclic!
assert_eq!(245.0, ds.calc_bias(0, 0)); assert_delta!(0.0, ds.calc_bias(0,0), 0.0000001);
} }
#[test] #[test]
fn calc_bias_offset_cyclic() { fn calc_biascyclic() {
let mut ds = build_hist_set(); let mut ds = build_hist_set();
ds.cyclic = true; ds.cyclic = true;
// 7th element -> x=7.5, x0=7.5 // 7th element -> x=7.5, x0=7.5
assert_eq!(0.0, ds.calc_bias(7, 0)); assert_delta!(1.0, ds.calc_bias(7, 0), 0.00000001);
// 8th element -> x=8.5, x0=7.5 // 8th element -> x=8.5, x0=7.5
assert_eq!(5.0, ds.calc_bias(8, 0)); assert_delta!(0.13472233779, ds.calc_bias(8, 0), 0.00000001);
// 1th element -> x=0.5, x0=7.5 // 1th element -> x=0.5, x0=7.5
// cyclic flag makes bin 0 neighboring bin 9, so the distance is actually 2 // cyclic flag makes bin 0 neighboring bin 9, so the distance is actually 2
assert_eq!(20.0, ds.calc_bias(0, 0)); assert_delta!(0.00032942643, ds.calc_bias(0, 0), 0.00000001);
// 2nd element -> x=1.5, x0=7.5 // 2nd element -> x=1.5, x0=7.5
assert_eq!(45.0, ds.calc_bias(1, 0)); assert_delta!(0.00000001, ds.calc_bias(1, 0), 0.00000001);
} }
#[test] #[test]

7
src/io.rs
View File

@@ -236,13 +236,6 @@ mod tests {
println!("{:?}", ds); println!("{:?}", ds);
assert_eq!(2, ds.num_windows); assert_eq!(2, ds.num_windows);
assert_eq!(cfg.num_bins[0], ds.dimens_lengths[0]); assert_eq!(cfg.num_bins[0], ds.dimens_lengths[0]);
// fields are private
// 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_pos);
// 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())
} }

318
src/lib.rs
View File

@@ -31,92 +31,52 @@ 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={:?} verbose={}, tolerance={}, iterations={}, temperature={}, cyclic={:?}", self.metadata_file, self.hist_min, self.hist_max, self.num_bins,
self.hist_max, self.num_bins, self.verbose, self.tolerance, self.verbose, self.tolerance, self.max_iterations, self.temperature, self.cyclic)
self.max_iterations, self.temperature, self.cyclic)
} }
} }
// 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 maximal difference for the calculated bias offsets is
// smaller then a tolerance value for every simulation window. // smaller then a tolerance value.
fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool { fn is_converged(old_F: &[f64], new_F: &[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 })
} }
// 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 given bias offsets (F)
// 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: &[f64]) -> f64 { fn calc_bin_probability(bin: usize, ds: &Dataset, F: &[f64]) -> f64 {
let mut denom_sum: f64 = 0.0; let mut denom_sum: f64 = 0.0;
let mut bin_count: f64 = 0.0; let mut bin_count: f64 = 0.0;
for (window, h) in ds.histograms.iter().enumerate() { for (window, h) in ds.histograms.iter().enumerate() {
bin_count += h.bins[bin]; bin_count += h.bins[bin];
let bias = ds.calc_bias(bin, window); let bias = ds.calc_bias(bin, window);
let bias_offset = ((F[window] - bias) / ds.kT).exp(); denom_sum += (h.num_points as f64) * bias * F[window];
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 and returns exp(F/kT)
fn calc_window_F(window: usize, ds: &Dataset, P: &[f64]) -> f64 { fn calc_window_F(window: usize, ds: &Dataset, P: &[f64]) -> f64 {
let bf_sum: f64 = (0..ds.num_bins).zip(P.iter()) // zip bins and P let f: f64 = (0..ds.num_bins).zip(P.iter()) // zip bins and P
.filter_map(|count_and_prob: (usize, &f64)| { .filter_map(|bin_and_prob: (usize, &f64)| {
if count_and_prob.1 == &0.0 { // skip zeros for speed if bin_and_prob.1 == &0.0 { // skip zeros for speed
None None
} else { } else {
Some(count_and_prob.1 * (-ds.calc_bias(count_and_prob.0, window) / ds.kT).exp()) let bias = ds.calc_bias(bin_and_prob.0, window);
Some(bin_and_prob.1 * bias)
} }
}).sum(); }).sum();
-ds.kT * bf_sum.ln() 1.0/f
} }
// 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: &[f64], F: &mut [f64], P: &mut [f64]) { fn perform_wham_iteration(ds: &Dataset, F_prev: &[f64], F: &mut [f64], P: &mut [f64]) {
// reset bias offsets
for window in 0..ds.num_windows {
F[window] = 0.0;
}
// for bin in 0..ds.num_bins {
// let x = get_x_for_bin(bin, ds.hist_min, ds.bin_width);
// let mut num = 0.0;
// let mut denom = 0.0;
// for window in 0..ds.num_windows {
// match ds.histograms[window].get_bin_count(bin) {
// Some(c) => num += c,
// _ => {}
// }
// let bias = calc_bias(
// ds.bias_fc[window],
// ds.bias_pos[window],
// x);
// let bf = ((F_prev[window]-bias) / ds.kT).exp();
// denom += ds.histograms[window].num_points as f64* bf
// }
// P[bin] = num / denom;
// for window in 0..ds.num_windows {
// let bias = calc_bias(
// ds.bias_fc[window],
// ds.bias_pos[window],
// x);
// let bf = (-bias/ds.kT).exp() * P[bin];
// F[window] += bf;
// }
// }
// for window in 0..ds.num_windows {
// F[window] = -ds.kT * F[window].ln();
// }
// evaluate first WHAM equation for each bin to // evaluate first WHAM equation for each bin to
// estimage probabilities based on previous offsets (F_prev) // estimage probabilities based on previous offsets (F_prev)
for bin in 0..ds.num_bins { for bin in 0..ds.num_bins {
@@ -128,9 +88,78 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &[f64], F: &mut [f64], P: &mut [
for window in 0..ds.num_windows { for window in 0..ds.num_windows {
F[window] = calc_window_F(window, ds, P); F[window] = calc_window_F(window, ds, P);
} }
} }
pub fn run(cfg: &Config) -> Result<(), Box<Error>>{
println!("Supplied WHAM options: {}", &cfg);
println!("Reading input files.");
// TODO Better error handling with nice error messages instead of a panic!
let histograms = io::read_data(&cfg)
.expect("No datapoints in histogram boundaries.");
println!("{}",&histograms);
// allocate required vectors.
let mut P: Vec<f64> = vec![f64::NAN; histograms.num_bins]; // bin probability
let mut F: Vec<f64> = vec![1.0; histograms.num_windows]; // bias offset exp(F/kT)
let mut F_prev: Vec<f64> = vec![f64::NAN; histograms.num_windows]; // previous bias offset
let mut F_tmp: Vec<f64> = vec![f64::NAN; histograms.num_windows]; // temp storage for F
let mut iteration = 0;
let mut converged = false;
// perform WHAM until convergence
while !converged && iteration < cfg.max_iterations {
iteration += 1;
// store F values before the next iteration
F_prev.copy_from_slice(&F);
// perform wham iteration (this updates F and P)
perform_wham_iteration(&histograms, &F_prev, &mut F, &mut P);
// convergence check
if iteration % 10 == 0 {
// This backups exp(F/kT) in a temporary vector and calculates true F and F_prev for
// convergence. Finally, F is restored. F_prev does not need to be restored because
// its overwritten for the next iteration.
F_tmp.copy_from_slice(&F);
for window in 0..histograms.num_windows {
F[window] = -histograms.kT * F[window].ln();
F_prev[window] = -histograms.kT * F_prev[window].ln();
}
converged = is_converged(&F_prev, &F, cfg.tolerance);
println!("Iteration {}: dF={}", &iteration, &diff_avg(&F_prev, &F));
F.copy_from_slice(&F_tmp);
}
// Dump free energy and bias offsets
//if iteration % 100 == 0 {
// free_energy(&histograms, &mut P, &mut A);
// dump_state(&histograms, &F, &F_prev, &P, &A);
//}
}
// Normalize P to sum(P) = 1.0
let P_sum: f64 = P.iter().sum();
P.iter_mut().map(|p| *p /= P_sum).count();
// calculate free energy and dump state
println!("Finished. Dumping final PMF");
let free_energy = calc_free_energy(&histograms, &P);
dump_state(&histograms, &F, &F_prev, &P, &free_energy);
if iteration == cfg.max_iterations {
println!("!!!!! WHAM not converged! (max iterations reached) !!!!!");
}
io::write_results(&cfg.output, &histograms, &free_energy, &P)?;
Ok(())
}
// get average difference between two bias offset sets // get average difference between two bias offset sets
fn diff_avg(F: &[f64], F_prev: &[f64]) -> f64 { fn diff_avg(F: &[f64], F_prev: &[f64]) -> f64 {
let mut F_sum: f64 = 0.0; let mut F_sum: f64 = 0.0;
@@ -138,90 +167,36 @@ fn diff_avg(F: &[f64], F_prev: &[f64]) -> f64 {
F_sum += (F[i]-F_prev[i]).abs() F_sum += (F[i]-F_prev[i]).abs()
} }
F_sum / F.len() as f64 F_sum / F.len() as f64
}
// calculate the normalized free energy from normalized probability values
fn free_energy(ds: &Dataset, P: &[f64], A: &mut [f64]) {
let mut bin_min = f64::MAX;
// Free energy calculation
for bin in 0..ds.num_bins {
A[bin] = -ds.kT*P[bin].ln();
if A[bin] < bin_min {
bin_min = A[bin];
}
}
// Make A relative to minimum
for bin in 0..ds.num_bins {
A[bin] -= bin_min;
}
} }
pub fn run(cfg: &Config) -> Result<(), Box<Error>>{ // calculate the normalized free energy from probability values
println!("Supplied WHAM options: {}", &cfg); fn calc_free_energy(ds: &Dataset, P: &[f64]) -> Vec<f64> {
let mut minimum = f64::MAX;
let mut free_energy: Vec<f64> = P.iter()
.map(|p| {
-ds.kT * p.ln()
})
.inspect(|free_e| {
if free_e < &minimum {
minimum = *free_e;
}
})
.collect();
// read input data into the histograms object for e in free_energy.iter_mut() {
println!("Reading input files."); *e -= minimum
}
let histograms = io::read_data(&cfg) // TODO nicer error handling for this free_energy
.expect("No datapoints in histogram boundaries.");
println!("{}",&histograms);
// allocate only once for better performance
let mut F_prev: Vec<f64> = vec![f64::INFINITY; histograms.num_windows];
let mut F: Vec<f64> = vec![0.0; histograms.num_windows];
let mut P: Vec<f64> = vec![f64::NAN; histograms.num_bins];
let mut A: Vec<f64> = vec![f64::NAN; histograms.num_bins];
// perform WHAM until convergence
let mut iteration = 0;
while !is_converged(&F_prev, &F, cfg.tolerance) && iteration < cfg.max_iterations {
iteration += 1;
// store F values before the next iteration
F_prev.copy_from_slice(&F[..]);
// perform wham iteration and update F
perform_wham_iteration(&histograms, &F_prev, &mut F, &mut P);
// output some stats during calculation
if iteration % 10 == 0 {
println!("Iteration {}: dF={}", &iteration, &diff_avg(&F_prev, &F));
}
// Dump free energy and bias offsets
if iteration % 100 == 0 {
free_energy(&histograms, &mut P, &mut A);
dump_state(&histograms, &F, &F_prev, &P, &A);
}
}
// Normalize P
let P_sum: f64 = P.iter().sum();
P.iter_mut().map(|p| *p /= P_sum).count();
// 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(())
} }
// TODO print nice headers for N dimensions
fn dump_state(ds: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], A: &[f64]) { fn dump_state(ds: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], A: &[f64]) {
let out = std::io::stdout(); let out = std::io::stdout();
let mut lock = out.lock(); let mut lock = out.lock();
writeln!(lock, "# PMF"); writeln!(lock, "# PMF");
writeln!(lock, "#x\t\tFree Energy\t\tP(x)"); writeln!(lock, "#x\t\tFree Energy\t\tP(x)");
for bin in 0..ds.num_bins { for bin in 0..ds.num_bins {
let x = ds.get_coords_for_bin(bin)[0]; // TODO let x = ds.get_coords_for_bin(bin)[0];
writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}", x, A[bin], P[bin]); writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}", x, A[bin], P[bin]);
} }
writeln!(lock, "# Bias offsets"); writeln!(lock, "# Bias offsets");
@@ -231,10 +206,26 @@ fn dump_state(ds: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], A: &[f64]) {
} }
} }
#[cfg(test)] #[cfg(test)]
mod tests { mod tests {
use super::histogram::{Dataset,Histogram}; use super::histogram::{Dataset,Histogram};
use std::f64; use std::f64;
use super::k_B;
macro_rules! assert_delta {
($x:expr, $y:expr, $d:expr) => {
assert!(($x-$y).abs() < $d, "{} != {}", $x, $y)
}
}
fn create_test_ds() -> Dataset {
let h1 = Histogram::new(10, vec![0.0, 1.0, 1.0, 8.0, 0.0]);
let h2 = Histogram::new(10, vec![0.0, 0.0, 8.0, 1.0, 1.0]);
Dataset::new(5, vec![5], vec![1.0], vec![0.0], vec![4.0],
vec![1.0, 1.0], vec![10.0, 10.0], 300.0*k_B, vec![h1, h2], false)
}
#[test] #[test]
fn is_converged() { fn is_converged() {
@@ -249,63 +240,44 @@ mod tests {
assert!(!converged); assert!(!converged);
} }
fn create_test_ds() -> Dataset { #[test]
let h1 = Histogram::new(10, vec![0.0, 0.0, 3.0, 4.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
let h2 = Histogram::new(20, vec![0.0, 0.0, 0.0, 3.0, 2.0, 5.0, 10.0, 0.0, 0.0, 0.0, 0.0]);
Dataset::new(4, vec![1], vec![1.0], vec![0.0], vec![4.0],
vec![1.0, 2.0], vec![10.0, 10.0], 2.479, vec![h1, h2], false)
}
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);
}
#[test]
fn calc_bias_offset() {
let ds = create_test_ds();
let probability = vec!(0.959, 0.331, 0.656, 46.750);
let expected = vec!(0.596, -0.250);
for window in 0..ds.num_windows {
let F = super::calc_window_F(window, &ds, &probability);
assert_near(expected[window], F, 0.001);
}
}
#[test]
#[ignore] // TODO
fn calc_bin_probability() { fn calc_bin_probability() {
let ds = create_test_ds(); let ds = create_test_ds();
let F = vec!(0.0, 0.0); let F = vec![1.0; ds.num_bins] ;
let expected = vec!(0.959, 0.331, 0.656, 46.750); let expected = vec!(0.0, 0.0825296687031316, 40.92355847097493,
for b in 0..4 { 124226.70003377, 2308526035.5283747);
for b in 0..ds.num_bins {
let p = super::calc_bin_probability(b, &ds, &F); let p = super::calc_bin_probability(b, &ds, &F);
assert_near(expected[b], p, 0.001); assert_delta!(expected[b], p, 0.0000001);
}
let F = vec!(1.0, 1.0);
let expected = vec!(0.641, 0.221, 0.439, 31.232);
for b in 0..4 {
let p = super::calc_bin_probability(b, &ds, &F);
assert_near(expected[b], p, 0.001);
} }
} }
#[test]
fn calc_bias_offset() {
let ds = create_test_ds();
let probability = vec!(0.0, 0.1, 0.2, 0.3, 0.4);
let expected = vec!(15.927477169990633, 15.927477169990633);
for window in 0..ds.num_windows {
let F = super::calc_window_F(window, &ds, &probability);
assert_delta!(expected[window], F, 0.0000001);
}
}
#[test] #[test]
#[ignore] // TODO
fn perform_wham_iteration() { fn perform_wham_iteration() {
let ds = create_test_ds(); let ds = create_test_ds();
let prev_F = vec![0.0; ds.num_windows]; let prev_F = vec![1.0; ds.num_windows];
let mut F = vec![0.0; ds.num_windows]; let mut F = vec![f64::NAN; ds.num_windows];
let mut P = vec![f64::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.5948, -0.2513); let expected_F = vec!(1.0, 1.0);
let expected_P = vec!(0.959, 0.331, 0.656, 46.750); let expected_P = vec!(0.0, 0.0825296687031316, 40.92355847097493,
124226.70003377, 2308526035.5283747);
for bin in 0..ds.num_bins { for bin in 0..ds.num_bins {
assert_near(expected_P[bin], P[bin], 0.01) assert_delta!(expected_P[bin], P[bin], 0.01)
} }
for window in 0..ds.num_windows { for window in 0..ds.num_windows {
assert_near(expected_F[window], F[window], 0.01) assert_delta!(expected_F[window], F[window], 0.01)
} }
} }