mirror of
https://github.com/dnlbauer/WHAM.git
synced 2026-09-10 22:25:31 +00:00
Compare commits
31 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
b6f338058a | ||
|
|
e0d2c1375a | ||
|
|
8802979a8e | ||
|
|
2d32d795f5 | ||
|
|
35ca1b1317 | ||
|
|
d9a219a36a | ||
|
|
2c565dee28 | ||
|
|
c05becea22 | ||
|
|
c27dbb7ab5 | ||
|
|
c9e18ebda9 | ||
|
|
0b5d0ebecd | ||
|
|
091b1f1382 | ||
|
|
9849c00321 | ||
|
|
034586d08e | ||
|
|
cc4448f1d8 | ||
|
|
5ca708abf1 | ||
|
|
31ffb5958d | ||
|
|
ffab55e0e1 | ||
|
|
69d0c7ca00 | ||
|
|
8a88d7742a | ||
|
18debf1a28 | ||
|
|
a856e5d34f | ||
|
|
2cc34919b0 | ||
|
|
cd57289211 | ||
|
|
9d8423bce9 | ||
|
|
2876c1dc3a | ||
|
|
b9d08ae24b | ||
|
|
07fc8344fe | ||
|
|
d7eea7aa03 | ||
|
|
3d93693eac | ||
|
|
38c30baa5b |
30
Cargo.lock
generated
30
Cargo.lock
generated
@@ -1,16 +1,5 @@
|
||||
# This file is automatically @generated by Cargo.
|
||||
# It is not intended for manual editing.
|
||||
[[package]]
|
||||
name = "GSL"
|
||||
version = "1.1.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "7830156ea389bcbbdc8f01bf140b609b892bf7cbd0ec6ccf9957ea2be6f25ad3"
|
||||
dependencies = [
|
||||
"c_vec",
|
||||
"libc",
|
||||
"pkg-config",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "addr2line"
|
||||
version = "0.13.0"
|
||||
@@ -78,12 +67,6 @@ version = "1.2.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "cf1de2fe8c75bc145a2f577add951f8134889b4795d47466a54a5c846d691693"
|
||||
|
||||
[[package]]
|
||||
name = "c_vec"
|
||||
version = "1.0.12"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "aa9e1d9f7d49e289f36f19effbf3d5a5e30163ecf9c7a3c9be94d5374dec5b9a"
|
||||
|
||||
[[package]]
|
||||
name = "cfg-if"
|
||||
version = "0.1.10"
|
||||
@@ -215,9 +198,9 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
|
||||
|
||||
[[package]]
|
||||
name = "libc"
|
||||
version = "0.2.79"
|
||||
version = "0.2.80"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "2448f6066e80e3bfc792e9c98bf705b4b0fc6e8ef5b43e5889aff0eaa9c58743"
|
||||
checksum = "4d58d1b70b004888f764dfbf6a26a3b0342a1632d33968e4a179d8011c760614"
|
||||
|
||||
[[package]]
|
||||
name = "memoffset"
|
||||
@@ -254,12 +237,6 @@ version = "0.21.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "37fd5004feb2ce328a52b0b3d01dbf4ffff72583493900ed15f22d4111c51693"
|
||||
|
||||
[[package]]
|
||||
name = "pkg-config"
|
||||
version = "0.3.19"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "3831453b3449ceb48b6d9c7ad7c96d5ea673e9b470a1dc578c2ce6521230884c"
|
||||
|
||||
[[package]]
|
||||
name = "ppv-lite86"
|
||||
version = "0.2.9"
|
||||
@@ -385,9 +362,8 @@ checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519"
|
||||
|
||||
[[package]]
|
||||
name = "wham"
|
||||
version = "0.9.9"
|
||||
version = "1.1.2"
|
||||
dependencies = [
|
||||
"GSL",
|
||||
"assert_approx_eq",
|
||||
"clap",
|
||||
"error-chain",
|
||||
|
||||
@@ -1,19 +1,21 @@
|
||||
[package]
|
||||
name = "wham"
|
||||
version = "0.9.9"
|
||||
version = "1.1.2"
|
||||
authors = ["Daniel Bauer <bauer@cbs.tu-darmstadt.de>"]
|
||||
description = "An implementation of the weighted histogram analysis method"
|
||||
license = "GPL-3.0"
|
||||
repository = "https://github.com/danijoo/WHAM"
|
||||
readme = "README.md"
|
||||
categories = ["science", "command-line-utilities", "molecular-dynamics", "algorithms"]
|
||||
categories = ["science", "command-line-utilities", "algorithms"]
|
||||
keywords = ["math", "statistics", "histogram", "bioinformatics", "molecular-dynamics"]
|
||||
exclude = [
|
||||
"example/*"
|
||||
]
|
||||
|
||||
[dependencies]
|
||||
clap = {version="2.32.0", features=['yaml']}
|
||||
error-chain = "0.12.0"
|
||||
rand = "0.7.*"
|
||||
GSL = "1.1"
|
||||
rayon = "1.0.3"
|
||||
|
||||
[dev-dependencies]
|
||||
|
||||
54
README.md
54
README.md
@@ -12,19 +12,14 @@ from umbrella sampling simulations. For more details on the method, I suggest *R
|
||||
Features
|
||||
---
|
||||
- Fast, especially for small systems
|
||||
- Multithreaded
|
||||
- Multidimensional
|
||||
- Error analysis
|
||||
- Multithreaded (automatically runs on all available cores)
|
||||
- Multidimensional (any number of collective variables are possible)
|
||||
- Autocorrelation to remove correlated samples
|
||||
- Error analysis via bootstrapping
|
||||
- Unit tested
|
||||
|
||||
Installation
|
||||
---
|
||||
WHAM requires the GSL library to be installed:
|
||||
```bash
|
||||
# on debian/ubuntu:
|
||||
sudo apt-get install libgsl0-dev
|
||||
```
|
||||
|
||||
Installation from source via cargo:
|
||||
```bash
|
||||
# cargo installation
|
||||
@@ -39,8 +34,8 @@ wham has a convenient command line interface. You can see all options with
|
||||
```wham -h```:
|
||||
|
||||
```
|
||||
wham 0.9.9
|
||||
D. Bauer <bauer@bio.tu-darmstadt.de>
|
||||
wham 1.1.0
|
||||
D. Bauer <bauer@cbs.tu-darmstadt.de>
|
||||
wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports
|
||||
potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
|
||||
|
||||
@@ -67,6 +62,8 @@ FLAGS:
|
||||
-c, --cyclic For periodic reaction coordinates. If this is set, the first and last coordinate bin in each
|
||||
dimension are treated as neighbors for the bias calculation.
|
||||
-h, --help Prints help information
|
||||
-g, --uncorr Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated
|
||||
samples (default is off).
|
||||
-V, --version Prints version information
|
||||
-v, --verbose Enables verbose output.
|
||||
|
||||
@@ -75,6 +72,10 @@ OPTIONS:
|
||||
--bt <bootstrap> Number of bayesian bootstrapping runs for error analysis by assigning random
|
||||
weights (defaults to 0).
|
||||
--seed <bootstrap_seed> Random seed for bootstrapping runs.
|
||||
--convdt <convdt> Performs WHAM for slices with the given delta in time and returns an output file
|
||||
for each slice. THis is useful to check the result for convergence. Example: with
|
||||
--convdt 100 and a timeseries ranging from 0-300, free energy surfaces for slices
|
||||
0-100, 0-200 and 0-300 will be given returned.
|
||||
--end <end> Skip rows in timeseries with an index larger than this value (defaults to 1e+20)
|
||||
-i, --iterations <ITERATIONS> Stop WHAM after this many iterations without convergence (defaults to 100,000).
|
||||
--max <HIST_MAX> Histogram maxima (comma separated). Also accepts "pi".
|
||||
@@ -85,10 +86,18 @@ OPTIONS:
|
||||
-T, --temperature <temperature> WHAM temperature in Kelvin.
|
||||
-t, --tolerance <TOLERANCE> Abortion criteria for WHAM calculation. WHAM stops if abs(F_new - F_old) <
|
||||
tolerance (defaults to 0.000001).
|
||||
|
||||
```
|
||||
|
||||
To run the two dimensional example (simulation of dialanine phi and psi angle):
|
||||
Examples
|
||||
---
|
||||
The example folder contains input and output files for two simple test systems:
|
||||
|
||||
- 1d_cyclic: Phi torsion angle of dialanine in vaccum
|
||||
- 2d_cyclic: Phi and psi torsion angles of the same system
|
||||
|
||||
The command below will run the two dimensional example (simulation of dialanine phi and psi angle) and calculate the free energy based on the two collective variables
|
||||
in the range of -3.14 to 3.14, with 100 bins in each dimension and periodic collective variables:
|
||||
|
||||
```bash
|
||||
wham --max 3.14,3.14 --min -3.14,-3.14 -T 300 --bins 100,100 --cyclic -f example/2d/metadata.dat
|
||||
> Supplied WHAM options: Metadata=example/2d/metadata.dat, hist_min=[-3.14, -3.14], hist_max=[3.14, 3.14], bins=[100, 100] verbose=false, tolerance=0.000001, iterations=100000, temperature=300, cyclic=true
|
||||
@@ -105,7 +114,6 @@ wham --max 3.14,3.14 --min -3.14,-3.14 -T 300 --bins 100,100 --cyclic -f example
|
||||
```
|
||||
After convergence, final bias offsets (F) and the free energy will be dumped to stdout and the output file is written.
|
||||
|
||||
|
||||
The output file contains the free energy and probability for each bin. Probabilities are normalized to sum to P=1.0 and
|
||||
the smallest free energy is set to 0 (with other free energies based on that).
|
||||
```
|
||||
@@ -125,7 +133,7 @@ the smallest free energy is set to 0 (with other free energies based on that).
|
||||
Error analysis
|
||||
---
|
||||
WHAM can perform error analysis using the bayesian bootstrapping method. Every simulation window is assumed to be an
|
||||
individual set of data point. By calculating probabilities N times with randomly assigned weights for each window,
|
||||
individual set of data points. By calculating probabilities N times with randomly assigned weights for each window,
|
||||
one can estimate the error as standard deviation between the N bootstrapping runs. For more details see
|
||||
*Van der Spoel, D. et al. (2010). g_wham—A Free Weighted Histogram Analysis Implementation Including Robust Error and
|
||||
Autocorrelation Estimates, JCTC, 6(12), 3713-3720*.
|
||||
@@ -134,17 +142,21 @@ To perform bayesian bootstrapping in WHAM, use the ```-bt <RUNS>``` flag to perf
|
||||
runs. The error estimates of bin probabilities and free energy will be given as standard error (SE) in a
|
||||
separate column (+/-) in the output file. If no error analysis is performed, these columns are set to 0.0.
|
||||
|
||||
Examples
|
||||
Autocorrelation analysis
|
||||
---
|
||||
The example folder contains input and output files for two simple test systems:
|
||||
With the ```--uncorr``` flag, WHAM calculates the autocorrelation time ```tau``` for all timeseries and all collective
|
||||
variables. Timeseries are then filtered based on their highest autocorrelation time to remove correlated samples from
|
||||
the dataset. This reduces the number of data points but can improve the accuracy of the result.
|
||||
|
||||
- 1d_cyclic: Phi torsion angle of dialanine in vaccum
|
||||
- 2d_cyclic: Phi and psi torsion angles of the same system
|
||||
For filtering, the statistical inefficiency `g` is calculated: ```g = 1 + 2*tau```, and only every `g`th element of the
|
||||
timeseries is used for unbiasing. A more detailed description of the method can be found in
|
||||
*Chodera, J.D. et al. (2007). Use of the weighted histogram analysis method for the analysis of simulated and parallel
|
||||
tempering simulations, JCTC 3(1):26-41*
|
||||
|
||||
|
||||
TODO
|
||||
---
|
||||
- Autocorrelation
|
||||
- Option to output histograms
|
||||
- Replica exchange
|
||||
|
||||
License & Citing
|
||||
@@ -152,7 +164,7 @@ License & Citing
|
||||
WHAM is licensed under the GPL-3.0 license. Please read the LICENSE file in this
|
||||
repository for more information.
|
||||
|
||||
There's no publication for this WHAM implementation. However, there is a citeabe DOI. If you use this software for your work, please consider citing it: *Bauer, D, WHAM - An efficient weighted histogram analysis implementation written in Rust, Zenodo. https://doi.org/10.5281/zenodo.1488597*
|
||||
There's no publication for this WHAM implementation. However, there is a citeabe DOI. If you use this software for your work, please consider citing it: *Bauer, D., WHAM - An efficient weighted histogram analysis implementation written in Rust, Zenodo. https://doi.org/10.5281/zenodo.1488597*
|
||||
|
||||
Parts of this work, especially some perfomance optimizations and the I/O format, are inspired by the
|
||||
implementation of A. Grossfield (*Grossfield, A, WHAM: the weighted histogram analysis method, http://membrane.urmc.rochester.edu/content/wham*).
|
||||
|
||||
@@ -1,101 +1,101 @@
|
||||
#coord1 Free Energy +/- Probability +/-
|
||||
-3.110177 7.158102 0.071438 0.003494 0.000068
|
||||
-3.047345 5.365727 0.069449 0.007168 0.000135
|
||||
-2.984513 3.873190 0.067495 0.013039 0.000232
|
||||
-2.921681 2.953162 0.067589 0.018855 0.000343
|
||||
-2.858849 1.949554 0.064640 0.028195 0.000480
|
||||
-2.796017 1.391747 0.063570 0.035261 0.000584
|
||||
-2.733186 1.128270 0.061710 0.039189 0.000620
|
||||
-2.670354 0.839970 0.060445 0.043991 0.000667
|
||||
-2.607522 0.624769 0.060762 0.047955 0.000739
|
||||
-2.544690 0.663757 0.060786 0.047211 0.000731
|
||||
-2.481858 1.051932 0.059774 0.040407 0.000617
|
||||
-2.419026 1.463048 0.060422 0.034267 0.000527
|
||||
-2.356194 1.990616 0.055242 0.027734 0.000360
|
||||
-2.293363 2.190692 0.044458 0.025597 0.000210
|
||||
-2.230531 2.553036 0.042144 0.022136 0.000159
|
||||
-2.167699 2.572522 0.043313 0.021964 0.000165
|
||||
-2.104867 2.472360 0.039799 0.022863 0.000157
|
||||
-2.042035 2.517562 0.036792 0.022453 0.000203
|
||||
-1.979203 2.469778 0.036324 0.022887 0.000237
|
||||
-1.916372 2.223125 0.049496 0.025266 0.000565
|
||||
-1.853540 2.080157 0.039205 0.026756 0.000453
|
||||
-1.790708 1.793841 0.041488 0.030011 0.000558
|
||||
-1.727876 1.458784 0.042362 0.034326 0.000565
|
||||
-1.665044 0.914441 0.037561 0.042697 0.000682
|
||||
-1.602212 0.303813 0.028701 0.054540 0.000761
|
||||
-1.539380 0.268335 0.025765 0.055321 0.000771
|
||||
-1.476549 0.000000 0.023371 0.061604 0.000866
|
||||
-1.413717 0.537061 0.026370 0.049671 0.000715
|
||||
-1.350885 1.439940 0.025620 0.034586 0.000547
|
||||
-1.288053 2.391838 0.027653 0.023614 0.000393
|
||||
-1.225221 3.779470 0.027093 0.013538 0.000229
|
||||
-1.162389 5.685555 0.027762 0.006305 0.000111
|
||||
-1.099557 7.661896 0.028996 0.002855 0.000051
|
||||
-1.036726 9.946594 0.032118 0.001142 0.000021
|
||||
-0.973894 12.359408 0.042801 0.000434 0.000009
|
||||
-0.911062 14.954655 0.052443 0.000153 0.000004
|
||||
-0.848230 17.742242 0.063261 0.000050 0.000002
|
||||
-0.785398 20.555785 0.068133 0.000016 0.000001
|
||||
-0.722566 22.811160 0.075518 0.000007 0.000000
|
||||
-0.659734 25.178094 0.086226 0.000003 0.000000
|
||||
-0.596903 26.442288 0.091783 0.000002 0.000000
|
||||
-0.534071 27.897565 0.093175 0.000001 0.000000
|
||||
-0.471239 29.062473 0.096027 0.000001 0.000000
|
||||
-0.408407 30.384521 0.096855 0.000000 0.000000
|
||||
-0.345575 31.638454 0.099956 0.000000 0.000000
|
||||
-0.282743 32.817727 0.107689 0.000000 0.000000
|
||||
-0.219911 33.770369 0.108946 0.000000 0.000000
|
||||
-0.157080 34.505503 0.112902 0.000000 0.000000
|
||||
-0.094248 35.431659 0.120856 0.000000 0.000000
|
||||
-0.031416 35.615810 0.122046 0.000000 0.000000
|
||||
0.031416 35.561946 0.118225 0.000000 0.000000
|
||||
0.094248 35.382089 0.108652 0.000000 0.000000
|
||||
0.157080 34.934827 0.118384 0.000000 0.000000
|
||||
0.219911 33.673460 0.119195 0.000000 0.000000
|
||||
0.282743 32.731563 0.121592 0.000000 0.000000
|
||||
0.345575 31.261855 0.129698 0.000000 0.000000
|
||||
0.408407 29.717377 0.139789 0.000000 0.000000
|
||||
0.471239 28.076620 0.137289 0.000001 0.000000
|
||||
0.534071 26.481097 0.135712 0.000002 0.000000
|
||||
0.596903 24.487358 0.135156 0.000003 0.000000
|
||||
0.659734 22.344251 0.131027 0.000008 0.000000
|
||||
0.722566 20.241543 0.133145 0.000018 0.000001
|
||||
0.785398 18.341869 0.131952 0.000039 0.000002
|
||||
0.848230 16.261582 0.134676 0.000091 0.000006
|
||||
0.911062 14.301801 0.134417 0.000199 0.000013
|
||||
0.973894 12.603788 0.131018 0.000394 0.000024
|
||||
1.036726 11.249601 0.131276 0.000678 0.000043
|
||||
1.099557 10.087886 0.132498 0.001079 0.000070
|
||||
1.162389 9.443303 0.132990 0.001398 0.000092
|
||||
1.225221 9.152799 0.132250 0.001570 0.000104
|
||||
1.288053 9.331937 0.133178 0.001462 0.000099
|
||||
1.350885 9.905546 0.133357 0.001161 0.000078
|
||||
1.413717 11.042050 0.133807 0.000736 0.000051
|
||||
1.476549 12.598167 0.132722 0.000395 0.000027
|
||||
1.539380 14.520167 0.131816 0.000183 0.000012
|
||||
1.602212 16.569783 0.131773 0.000080 0.000005
|
||||
1.665044 18.687390 0.132601 0.000034 0.000002
|
||||
1.727876 20.775408 0.133867 0.000015 0.000001
|
||||
1.790708 22.905200 0.129265 0.000006 0.000000
|
||||
1.853540 24.643852 0.128894 0.000003 0.000000
|
||||
1.916372 26.301740 0.130266 0.000002 0.000000
|
||||
1.979203 27.372071 0.128620 0.000001 0.000000
|
||||
2.042035 28.697726 0.133263 0.000001 0.000000
|
||||
2.104867 29.417901 0.133513 0.000000 0.000000
|
||||
2.167699 30.008351 0.130925 0.000000 0.000000
|
||||
2.230531 30.406016 0.124398 0.000000 0.000000
|
||||
2.293363 30.171275 0.122493 0.000000 0.000000
|
||||
2.356194 29.884646 0.130278 0.000000 0.000000
|
||||
2.419026 29.428153 0.130844 0.000000 0.000000
|
||||
2.481858 28.546982 0.148114 0.000001 0.000000
|
||||
2.544690 27.757520 0.133045 0.000001 0.000000
|
||||
2.607522 26.505787 0.134364 0.000001 0.000000
|
||||
2.670354 24.491866 0.115061 0.000003 0.000000
|
||||
2.733186 22.320664 0.110036 0.000008 0.000000
|
||||
2.796017 20.052723 0.107894 0.000020 0.000001
|
||||
2.858849 17.655650 0.105964 0.000052 0.000002
|
||||
2.921681 15.471590 0.107005 0.000125 0.000005
|
||||
2.984513 13.138167 0.099378 0.000318 0.000012
|
||||
3.047345 11.092386 0.087770 0.000722 0.000022
|
||||
3.110177 9.065722 0.077092 0.001626 0.000038
|
||||
-3.110177 7.158102 0.000000 0.003494 0.000000
|
||||
-3.047345 5.365727 0.000000 0.007168 0.000000
|
||||
-2.984513 3.873190 0.000000 0.013039 0.000000
|
||||
-2.921681 2.953162 0.000000 0.018855 0.000000
|
||||
-2.858849 1.949554 0.000000 0.028195 0.000000
|
||||
-2.796017 1.391747 0.000000 0.035261 0.000000
|
||||
-2.733186 1.128270 0.000000 0.039189 0.000000
|
||||
-2.670354 0.839970 0.000000 0.043991 0.000000
|
||||
-2.607522 0.624769 0.000000 0.047955 0.000000
|
||||
-2.544690 0.663757 0.000000 0.047211 0.000000
|
||||
-2.481858 1.051932 0.000000 0.040407 0.000000
|
||||
-2.419026 1.463048 0.000000 0.034267 0.000000
|
||||
-2.356194 1.990616 0.000000 0.027734 0.000000
|
||||
-2.293363 2.190692 0.000000 0.025597 0.000000
|
||||
-2.230531 2.553036 0.000000 0.022136 0.000000
|
||||
-2.167699 2.572522 0.000000 0.021964 0.000000
|
||||
-2.104867 2.472360 0.000000 0.022863 0.000000
|
||||
-2.042035 2.517562 0.000000 0.022453 0.000000
|
||||
-1.979203 2.469778 0.000000 0.022887 0.000000
|
||||
-1.916372 2.223125 0.000000 0.025266 0.000000
|
||||
-1.853540 2.080157 0.000000 0.026756 0.000000
|
||||
-1.790708 1.793841 0.000000 0.030011 0.000000
|
||||
-1.727876 1.458784 0.000000 0.034326 0.000000
|
||||
-1.665044 0.914441 0.000000 0.042697 0.000000
|
||||
-1.602212 0.303813 0.000000 0.054540 0.000000
|
||||
-1.539380 0.268335 0.000000 0.055321 0.000000
|
||||
-1.476549 0.000000 0.000000 0.061604 0.000000
|
||||
-1.413717 0.537061 0.000000 0.049671 0.000000
|
||||
-1.350885 1.439940 0.000000 0.034586 0.000000
|
||||
-1.288053 2.391838 0.000000 0.023614 0.000000
|
||||
-1.225221 3.779470 0.000000 0.013538 0.000000
|
||||
-1.162389 5.685555 0.000000 0.006305 0.000000
|
||||
-1.099557 7.661896 0.000000 0.002855 0.000000
|
||||
-1.036726 9.946594 0.000000 0.001142 0.000000
|
||||
-0.973894 12.359408 0.000000 0.000434 0.000000
|
||||
-0.911062 14.954655 0.000000 0.000153 0.000000
|
||||
-0.848230 17.742242 0.000000 0.000050 0.000000
|
||||
-0.785398 20.555785 0.000000 0.000016 0.000000
|
||||
-0.722566 22.811160 0.000000 0.000007 0.000000
|
||||
-0.659734 25.178094 0.000000 0.000003 0.000000
|
||||
-0.596903 26.442288 0.000000 0.000002 0.000000
|
||||
-0.534071 27.897565 0.000000 0.000001 0.000000
|
||||
-0.471239 29.062473 0.000000 0.000001 0.000000
|
||||
-0.408407 30.384521 0.000000 0.000000 0.000000
|
||||
-0.345575 31.638454 0.000000 0.000000 0.000000
|
||||
-0.282743 32.817727 0.000000 0.000000 0.000000
|
||||
-0.219911 33.770369 0.000000 0.000000 0.000000
|
||||
-0.157080 34.505503 0.000000 0.000000 0.000000
|
||||
-0.094248 35.431659 0.000000 0.000000 0.000000
|
||||
-0.031416 35.615810 0.000000 0.000000 0.000000
|
||||
0.031416 35.561946 0.000000 0.000000 0.000000
|
||||
0.094248 35.382089 0.000000 0.000000 0.000000
|
||||
0.157080 34.934827 0.000000 0.000000 0.000000
|
||||
0.219911 33.673460 0.000000 0.000000 0.000000
|
||||
0.282743 32.731563 0.000000 0.000000 0.000000
|
||||
0.345575 31.261855 0.000000 0.000000 0.000000
|
||||
0.408407 29.717377 0.000000 0.000000 0.000000
|
||||
0.471239 28.076620 0.000000 0.000001 0.000000
|
||||
0.534071 26.481097 0.000000 0.000002 0.000000
|
||||
0.596903 24.487358 0.000000 0.000003 0.000000
|
||||
0.659734 22.344251 0.000000 0.000008 0.000000
|
||||
0.722566 20.241543 0.000000 0.000018 0.000000
|
||||
0.785398 18.341869 0.000000 0.000039 0.000000
|
||||
0.848230 16.261582 0.000000 0.000091 0.000000
|
||||
0.911062 14.301801 0.000000 0.000199 0.000000
|
||||
0.973894 12.603788 0.000000 0.000394 0.000000
|
||||
1.036726 11.249601 0.000000 0.000678 0.000000
|
||||
1.099557 10.087886 0.000000 0.001079 0.000000
|
||||
1.162389 9.443303 0.000000 0.001398 0.000000
|
||||
1.225221 9.152799 0.000000 0.001570 0.000000
|
||||
1.288053 9.331937 0.000000 0.001462 0.000000
|
||||
1.350885 9.905546 0.000000 0.001161 0.000000
|
||||
1.413717 11.042050 0.000000 0.000736 0.000000
|
||||
1.476549 12.598167 0.000000 0.000395 0.000000
|
||||
1.539380 14.520167 0.000000 0.000183 0.000000
|
||||
1.602212 16.569783 0.000000 0.000080 0.000000
|
||||
1.665044 18.687390 0.000000 0.000034 0.000000
|
||||
1.727876 20.775408 0.000000 0.000015 0.000000
|
||||
1.790708 22.905200 0.000000 0.000006 0.000000
|
||||
1.853540 24.643852 0.000000 0.000003 0.000000
|
||||
1.916372 26.301740 0.000000 0.000002 0.000000
|
||||
1.979203 27.372071 0.000000 0.000001 0.000000
|
||||
2.042035 28.697726 0.000000 0.000001 0.000000
|
||||
2.104867 29.417901 0.000000 0.000000 0.000000
|
||||
2.167699 30.008351 0.000000 0.000000 0.000000
|
||||
2.230531 30.406016 0.000000 0.000000 0.000000
|
||||
2.293363 30.171275 0.000000 0.000000 0.000000
|
||||
2.356194 29.884646 0.000000 0.000000 0.000000
|
||||
2.419026 29.428153 0.000000 0.000000 0.000000
|
||||
2.481858 28.546982 0.000000 0.000001 0.000000
|
||||
2.544690 27.757520 0.000000 0.000001 0.000000
|
||||
2.607522 26.505787 0.000000 0.000001 0.000000
|
||||
2.670354 24.491866 0.000000 0.000003 0.000000
|
||||
2.733186 22.320664 0.000000 0.000008 0.000000
|
||||
2.796017 20.052723 0.000000 0.000020 0.000000
|
||||
2.858849 17.655650 0.000000 0.000052 0.000000
|
||||
2.921681 15.471590 0.000000 0.000125 0.000000
|
||||
2.984513 13.138167 0.000000 0.000318 0.000000
|
||||
3.047345 11.092386 0.000000 0.000722 0.000000
|
||||
3.110177 9.065722 0.000000 0.001626 0.000000
|
||||
|
||||
101
example/1d_cyclic/wham_bt.out
Normal file
101
example/1d_cyclic/wham_bt.out
Normal file
@@ -0,0 +1,101 @@
|
||||
#coord1 Free Energy +/- Probability +/-
|
||||
-3.110177 7.158102 0.071438 0.003494 0.000068
|
||||
-3.047345 5.365727 0.069449 0.007168 0.000135
|
||||
-2.984513 3.873190 0.067495 0.013039 0.000232
|
||||
-2.921681 2.953162 0.067589 0.018855 0.000343
|
||||
-2.858849 1.949554 0.064640 0.028195 0.000480
|
||||
-2.796017 1.391747 0.063570 0.035261 0.000584
|
||||
-2.733186 1.128270 0.061710 0.039189 0.000620
|
||||
-2.670354 0.839970 0.060445 0.043991 0.000667
|
||||
-2.607522 0.624769 0.060762 0.047955 0.000739
|
||||
-2.544690 0.663757 0.060786 0.047211 0.000731
|
||||
-2.481858 1.051932 0.059774 0.040407 0.000617
|
||||
-2.419026 1.463048 0.060422 0.034267 0.000527
|
||||
-2.356194 1.990616 0.055242 0.027734 0.000360
|
||||
-2.293363 2.190692 0.044458 0.025597 0.000210
|
||||
-2.230531 2.553036 0.042144 0.022136 0.000159
|
||||
-2.167699 2.572522 0.043313 0.021964 0.000165
|
||||
-2.104867 2.472360 0.039799 0.022863 0.000157
|
||||
-2.042035 2.517562 0.036792 0.022453 0.000203
|
||||
-1.979203 2.469778 0.036324 0.022887 0.000237
|
||||
-1.916372 2.223125 0.049496 0.025266 0.000565
|
||||
-1.853540 2.080157 0.039205 0.026756 0.000453
|
||||
-1.790708 1.793841 0.041488 0.030011 0.000558
|
||||
-1.727876 1.458784 0.042362 0.034326 0.000565
|
||||
-1.665044 0.914441 0.037561 0.042697 0.000682
|
||||
-1.602212 0.303813 0.028701 0.054540 0.000761
|
||||
-1.539380 0.268335 0.025765 0.055321 0.000771
|
||||
-1.476549 0.000000 0.023371 0.061604 0.000866
|
||||
-1.413717 0.537061 0.026370 0.049671 0.000715
|
||||
-1.350885 1.439940 0.025620 0.034586 0.000547
|
||||
-1.288053 2.391838 0.027653 0.023614 0.000393
|
||||
-1.225221 3.779470 0.027093 0.013538 0.000229
|
||||
-1.162389 5.685555 0.027762 0.006305 0.000111
|
||||
-1.099557 7.661896 0.028996 0.002855 0.000051
|
||||
-1.036726 9.946594 0.032118 0.001142 0.000021
|
||||
-0.973894 12.359408 0.042801 0.000434 0.000009
|
||||
-0.911062 14.954655 0.052443 0.000153 0.000004
|
||||
-0.848230 17.742242 0.063261 0.000050 0.000002
|
||||
-0.785398 20.555785 0.068133 0.000016 0.000001
|
||||
-0.722566 22.811160 0.075518 0.000007 0.000000
|
||||
-0.659734 25.178094 0.086226 0.000003 0.000000
|
||||
-0.596903 26.442288 0.091783 0.000002 0.000000
|
||||
-0.534071 27.897565 0.093175 0.000001 0.000000
|
||||
-0.471239 29.062473 0.096027 0.000001 0.000000
|
||||
-0.408407 30.384521 0.096855 0.000000 0.000000
|
||||
-0.345575 31.638454 0.099956 0.000000 0.000000
|
||||
-0.282743 32.817727 0.107689 0.000000 0.000000
|
||||
-0.219911 33.770369 0.108946 0.000000 0.000000
|
||||
-0.157080 34.505503 0.112902 0.000000 0.000000
|
||||
-0.094248 35.431659 0.120856 0.000000 0.000000
|
||||
-0.031416 35.615810 0.122046 0.000000 0.000000
|
||||
0.031416 35.561946 0.118225 0.000000 0.000000
|
||||
0.094248 35.382089 0.108652 0.000000 0.000000
|
||||
0.157080 34.934827 0.118384 0.000000 0.000000
|
||||
0.219911 33.673460 0.119195 0.000000 0.000000
|
||||
0.282743 32.731563 0.121592 0.000000 0.000000
|
||||
0.345575 31.261855 0.129698 0.000000 0.000000
|
||||
0.408407 29.717377 0.139789 0.000000 0.000000
|
||||
0.471239 28.076620 0.137289 0.000001 0.000000
|
||||
0.534071 26.481097 0.135712 0.000002 0.000000
|
||||
0.596903 24.487358 0.135156 0.000003 0.000000
|
||||
0.659734 22.344251 0.131027 0.000008 0.000000
|
||||
0.722566 20.241543 0.133145 0.000018 0.000001
|
||||
0.785398 18.341869 0.131952 0.000039 0.000002
|
||||
0.848230 16.261582 0.134676 0.000091 0.000006
|
||||
0.911062 14.301801 0.134417 0.000199 0.000013
|
||||
0.973894 12.603788 0.131018 0.000394 0.000024
|
||||
1.036726 11.249601 0.131276 0.000678 0.000043
|
||||
1.099557 10.087886 0.132498 0.001079 0.000070
|
||||
1.162389 9.443303 0.132990 0.001398 0.000092
|
||||
1.225221 9.152799 0.132250 0.001570 0.000104
|
||||
1.288053 9.331937 0.133178 0.001462 0.000099
|
||||
1.350885 9.905546 0.133357 0.001161 0.000078
|
||||
1.413717 11.042050 0.133807 0.000736 0.000051
|
||||
1.476549 12.598167 0.132722 0.000395 0.000027
|
||||
1.539380 14.520167 0.131816 0.000183 0.000012
|
||||
1.602212 16.569783 0.131773 0.000080 0.000005
|
||||
1.665044 18.687390 0.132601 0.000034 0.000002
|
||||
1.727876 20.775408 0.133867 0.000015 0.000001
|
||||
1.790708 22.905200 0.129265 0.000006 0.000000
|
||||
1.853540 24.643852 0.128894 0.000003 0.000000
|
||||
1.916372 26.301740 0.130266 0.000002 0.000000
|
||||
1.979203 27.372071 0.128620 0.000001 0.000000
|
||||
2.042035 28.697726 0.133263 0.000001 0.000000
|
||||
2.104867 29.417901 0.133513 0.000000 0.000000
|
||||
2.167699 30.008351 0.130925 0.000000 0.000000
|
||||
2.230531 30.406016 0.124398 0.000000 0.000000
|
||||
2.293363 30.171275 0.122493 0.000000 0.000000
|
||||
2.356194 29.884646 0.130278 0.000000 0.000000
|
||||
2.419026 29.428153 0.130844 0.000000 0.000000
|
||||
2.481858 28.546982 0.148114 0.000001 0.000000
|
||||
2.544690 27.757520 0.133045 0.000001 0.000000
|
||||
2.607522 26.505787 0.134364 0.000001 0.000000
|
||||
2.670354 24.491866 0.115061 0.000003 0.000000
|
||||
2.733186 22.320664 0.110036 0.000008 0.000000
|
||||
2.796017 20.052723 0.107894 0.000020 0.000001
|
||||
2.858849 17.655650 0.105964 0.000052 0.000002
|
||||
2.921681 15.471590 0.107005 0.000125 0.000005
|
||||
2.984513 13.138167 0.099378 0.000318 0.000012
|
||||
3.047345 11.092386 0.087770 0.000722 0.000022
|
||||
3.110177 9.065722 0.077092 0.001626 0.000038
|
||||
101
example/1d_cyclic/wham_uncorrelated.out
Normal file
101
example/1d_cyclic/wham_uncorrelated.out
Normal file
@@ -0,0 +1,101 @@
|
||||
#coord1 Free Energy +/- Probability +/-
|
||||
-3.110177 7.531315 0.000000 0.003080 0.000000
|
||||
-3.047345 5.690157 0.000000 0.006443 0.000000
|
||||
-2.984513 4.243063 0.000000 0.011509 0.000000
|
||||
-2.921681 3.334686 0.000000 0.016564 0.000000
|
||||
-2.858849 2.277349 0.000000 0.025309 0.000000
|
||||
-2.796017 1.723296 0.000000 0.031604 0.000000
|
||||
-2.733186 1.246264 0.000000 0.038265 0.000000
|
||||
-2.670354 1.099867 0.000000 0.040578 0.000000
|
||||
-2.607522 0.771910 0.000000 0.046279 0.000000
|
||||
-2.544690 0.770616 0.000000 0.046303 0.000000
|
||||
-2.481858 1.265507 0.000000 0.037971 0.000000
|
||||
-2.419026 1.562335 0.000000 0.033711 0.000000
|
||||
-2.356194 1.891577 0.000000 0.029542 0.000000
|
||||
-2.293363 2.227858 0.000000 0.025816 0.000000
|
||||
-2.230531 2.488355 0.000000 0.023256 0.000000
|
||||
-2.167699 2.502265 0.000000 0.023127 0.000000
|
||||
-2.104867 2.358037 0.000000 0.024503 0.000000
|
||||
-2.042035 2.278147 0.000000 0.025301 0.000000
|
||||
-1.979203 2.974067 0.000000 0.019141 0.000000
|
||||
-1.916372 2.696600 0.000000 0.021393 0.000000
|
||||
-1.853540 2.361827 0.000000 0.024466 0.000000
|
||||
-1.790708 1.516746 0.000000 0.034332 0.000000
|
||||
-1.727876 1.526829 0.000000 0.034194 0.000000
|
||||
-1.665044 0.884114 0.000000 0.044244 0.000000
|
||||
-1.602212 0.323912 0.000000 0.055385 0.000000
|
||||
-1.539380 0.197985 0.000000 0.058253 0.000000
|
||||
-1.476549 0.000000 0.000000 0.063065 0.000000
|
||||
-1.413717 0.458247 0.000000 0.052481 0.000000
|
||||
-1.350885 1.389410 0.000000 0.036131 0.000000
|
||||
-1.288053 2.386522 0.000000 0.024225 0.000000
|
||||
-1.225221 3.743253 0.000000 0.014062 0.000000
|
||||
-1.162389 5.566654 0.000000 0.006770 0.000000
|
||||
-1.099557 7.822800 0.000000 0.002740 0.000000
|
||||
-1.036726 10.128719 0.000000 0.001087 0.000000
|
||||
-0.973894 12.199246 0.000000 0.000474 0.000000
|
||||
-0.911062 14.488129 0.000000 0.000189 0.000000
|
||||
-0.848230 16.902310 0.000000 0.000072 0.000000
|
||||
-0.785398 18.910200 0.000000 0.000032 0.000000
|
||||
-0.722566 21.241681 0.000000 0.000013 0.000000
|
||||
-0.659734 22.706373 0.000000 0.000007 0.000000
|
||||
-0.596903 24.531129 0.000000 0.000003 0.000000
|
||||
-0.534071 25.936227 0.000000 0.000002 0.000000
|
||||
-0.471239 27.000262 0.000000 0.000001 0.000000
|
||||
-0.408407 28.673293 0.000000 0.000001 0.000000
|
||||
-0.345575 29.335203 0.000000 0.000000 0.000000
|
||||
-0.282743 30.841118 0.000000 0.000000 0.000000
|
||||
-0.219911 31.983859 0.000000 0.000000 0.000000
|
||||
-0.157080 32.144015 0.000000 0.000000 0.000000
|
||||
-0.094248 33.885395 0.000000 0.000000 0.000000
|
||||
-0.031416 33.783105 0.000000 0.000000 0.000000
|
||||
0.031416 34.243727 0.000000 0.000000 0.000000
|
||||
0.094248 33.975567 0.000000 0.000000 0.000000
|
||||
0.157080 32.994787 0.000000 0.000000 0.000000
|
||||
0.219911 32.607398 0.000000 0.000000 0.000000
|
||||
0.282743 31.401902 0.000000 0.000000 0.000000
|
||||
0.345575 29.911670 0.000000 0.000000 0.000000
|
||||
0.408407 28.603574 0.000000 0.000001 0.000000
|
||||
0.471239 26.925443 0.000000 0.000001 0.000000
|
||||
0.534071 25.298070 0.000000 0.000002 0.000000
|
||||
0.596903 23.638560 0.000000 0.000005 0.000000
|
||||
0.659734 21.156231 0.000000 0.000013 0.000000
|
||||
0.722566 19.126480 0.000000 0.000029 0.000000
|
||||
0.785398 17.351953 0.000000 0.000060 0.000000
|
||||
0.848230 15.135525 0.000000 0.000146 0.000000
|
||||
0.911062 13.188112 0.000000 0.000319 0.000000
|
||||
0.973894 11.536983 0.000000 0.000618 0.000000
|
||||
1.036726 10.158328 0.000000 0.001074 0.000000
|
||||
1.099557 9.109036 0.000000 0.001636 0.000000
|
||||
1.162389 8.282343 0.000000 0.002279 0.000000
|
||||
1.225221 8.022102 0.000000 0.002530 0.000000
|
||||
1.288053 8.162415 0.000000 0.002391 0.000000
|
||||
1.350885 8.600135 0.000000 0.002006 0.000000
|
||||
1.413717 9.837348 0.000000 0.001222 0.000000
|
||||
1.476549 11.363156 0.000000 0.000663 0.000000
|
||||
1.539380 13.077849 0.000000 0.000333 0.000000
|
||||
1.602212 15.353594 0.000000 0.000134 0.000000
|
||||
1.665044 17.565051 0.000000 0.000055 0.000000
|
||||
1.727876 19.710884 0.000000 0.000023 0.000000
|
||||
1.790708 21.721260 0.000000 0.000010 0.000000
|
||||
1.853540 23.567649 0.000000 0.000005 0.000000
|
||||
1.916372 25.008817 0.000000 0.000003 0.000000
|
||||
1.979203 26.405367 0.000000 0.000002 0.000000
|
||||
2.042035 28.070821 0.000000 0.000001 0.000000
|
||||
2.104867 28.877213 0.000000 0.000001 0.000000
|
||||
2.167699 29.378146 0.000000 0.000000 0.000000
|
||||
2.230531 31.093267 0.000000 0.000000 0.000000
|
||||
2.293363 30.704994 0.000000 0.000000 0.000000
|
||||
2.356194 30.563093 0.000000 0.000000 0.000000
|
||||
2.419026 31.215952 0.000000 0.000000 0.000000
|
||||
2.481858 30.331416 0.000000 0.000000 0.000000
|
||||
2.544690 29.005123 0.000000 0.000001 0.000000
|
||||
2.607522 27.674618 0.000000 0.000001 0.000000
|
||||
2.670354 24.911788 0.000000 0.000003 0.000000
|
||||
2.733186 22.746417 0.000000 0.000007 0.000000
|
||||
2.796017 20.608288 0.000000 0.000016 0.000000
|
||||
2.858849 18.018280 0.000000 0.000046 0.000000
|
||||
2.921681 15.949215 0.000000 0.000105 0.000000
|
||||
2.984513 13.617809 0.000000 0.000268 0.000000
|
||||
3.047345 11.432193 0.000000 0.000645 0.000000
|
||||
3.110177 9.461494 0.000000 0.001420 0.000000
|
||||
20
src/cli.yml
20
src/cli.yml
@@ -1,6 +1,6 @@
|
||||
name: wham
|
||||
version: "0.9.9"
|
||||
author: D. Bauer <bauer@bio.tu-darmstadt.de>
|
||||
version: "1.1.2"
|
||||
author: D. Bauer <bauer@cbs.tu-darmstadt.de>
|
||||
about: |
|
||||
wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
|
||||
|
||||
@@ -101,3 +101,19 @@ args:
|
||||
help: Skip rows in timeseries with an index larger than this value (defaults to 1e+20)
|
||||
takes_value: true
|
||||
required: false
|
||||
- uncorr:
|
||||
short: g
|
||||
long: uncorr
|
||||
help: Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated samples (default is off).
|
||||
takes_value: false
|
||||
required: false
|
||||
- convdt:
|
||||
long: convdt
|
||||
help: "Performs WHAM for slices with the given delta in time and returns an output file for each slice. THis is useful to check the result for convergence. Example: with --convdt 100 and a timeseries ranging from 0-300, free energy surfaces for slices 0-100, 0-200 and 0-300 will be given returned."
|
||||
takes_value: true
|
||||
required: false
|
||||
- ignore_empty:
|
||||
long: ignore_empty
|
||||
help: If this is set, do not fail if a histogram is empty.
|
||||
takes_value: false
|
||||
required: false
|
||||
|
||||
78
src/correlation_analysis.rs
Normal file
78
src/correlation_analysis.rs
Normal file
@@ -0,0 +1,78 @@
|
||||
use super::statistics;
|
||||
|
||||
// calculates the statistical inefficiency g of the given timeseries
|
||||
// the quantity g can be thought of: N/g is the number of uncorrelated
|
||||
// configurations in the timeseries, where samples are separated by
|
||||
// the a multiple of g
|
||||
// For details, see "Chodera et al. (2007). Use of a Weighted Histogram Analysis
|
||||
// Method for the Analysis of Simulated and Parallel Tempering Simulations, JCTC"
|
||||
pub fn statistical_ineff(timeseries: &[f64]) -> f64 {
|
||||
let n = timeseries.len();
|
||||
let mean = statistics::mean(timeseries);
|
||||
let d_mean = timeseries.iter().map(|x| x-mean).collect::<Vec<f64>>();
|
||||
let cov = statistics::autocov(timeseries);
|
||||
|
||||
let mut g = 1.0;
|
||||
for t in 1..(n-1) {
|
||||
let tmp = d_mean[0..n-t].iter().zip(d_mean[t..n].iter()).map(|(x,y)| x*y);
|
||||
let c = tmp.map(|x| x+x).sum::<f64>() / (2.0 * (n as f64 - t as f64)*cov);
|
||||
if c <= 0.0 {
|
||||
break;
|
||||
}
|
||||
g = g + (2.0*c*(1.0-t as f64/n as f64))
|
||||
}
|
||||
if g < 1.0 {
|
||||
1.0
|
||||
} else {
|
||||
g
|
||||
}
|
||||
}
|
||||
|
||||
// The autocorrelation time of a timeseries can be deduced from the
|
||||
// `statistical_ineff` by (g-1)/2.0
|
||||
pub fn autocorrelation_time(g: f64) -> f64 {
|
||||
(g - 1.0) / 2.0
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use std::io::{BufRead, BufReader};
|
||||
use std::fs::File;
|
||||
|
||||
fn read_timeseries(filename: &str) -> Vec<f64> {
|
||||
let mut timeseries: Vec<f64> = Vec::new();
|
||||
let file = File::open(filename).unwrap();
|
||||
let reader = BufReader::new(&file);
|
||||
for line in reader.lines() {
|
||||
let val = line.unwrap().split_whitespace()
|
||||
.collect::<Vec<&str>>()[1].parse::<f64>().unwrap();
|
||||
timeseries.push(val);
|
||||
}
|
||||
timeseries
|
||||
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn statistical_ineff() {
|
||||
let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
|
||||
let g = super::statistical_ineff(×eries);
|
||||
println!("{:?}", g);
|
||||
assert!((g - 3.859).abs() < 0.001);
|
||||
|
||||
// a "random" timeseries with g < 1.0
|
||||
let timeseries = [1_f64, 4_f64, 921_f64, 121213_f64, 23192_f64,
|
||||
8913_f64, 1232_f64, 2_f64, 151_f64, 123091_f64];
|
||||
let g = super::statistical_ineff(×eries);
|
||||
assert_approx_eq!(g, 1.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn autocorrelation_time() {
|
||||
let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
|
||||
let g = super::statistical_ineff(×eries);
|
||||
let tau = super::autocorrelation_time(g);
|
||||
println!("{:?}", tau);
|
||||
assert!((tau - 1.430).abs() < 0.001)
|
||||
}
|
||||
|
||||
}
|
||||
@@ -2,7 +2,7 @@ use rand::prelude::*;
|
||||
use super::histogram::{Dataset};
|
||||
use super::perform_wham;
|
||||
use super::{Config,calc_free_energy};
|
||||
use rgsl::statistics;
|
||||
use super::statistics;
|
||||
|
||||
// returns a set of num_windows continious weights by
|
||||
// a) generate num_windows-1 random variables and sort them
|
||||
@@ -48,7 +48,7 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec<f64>,Ve
|
||||
let mut P_se = vec![0.0; ds.num_bins];
|
||||
for bin in 0..ds.num_bins {
|
||||
let Ps = bootstrapped_Ps.iter().map(|window| window[bin]).collect::<Vec<f64>>();
|
||||
P_se[bin] = statistics::sd(&Ps, 1, num_runs)/(num_runs as f64).sqrt();
|
||||
P_se[bin] = statistics::sd(&Ps)/(num_runs as f64).sqrt();
|
||||
}
|
||||
|
||||
// SE of A
|
||||
@@ -60,13 +60,13 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec<f64>,Ve
|
||||
let mut A_se = vec![0.0; ds.num_bins];
|
||||
for bin in 0..ds.num_bins {
|
||||
let As = bootstrapped_As.iter().map(|window| window[bin]).collect::<Vec<f64>>();
|
||||
A_se[bin] = statistics::sd(&As, 1, num_runs)/(num_runs as f64).sqrt();
|
||||
A_se[bin] = statistics::sd(&As)/(num_runs as f64).sqrt();
|
||||
}
|
||||
|
||||
(P_se, A_se)
|
||||
}
|
||||
|
||||
#[cfg(tests)]
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use super::super::k_B;
|
||||
@@ -99,8 +99,9 @@ mod tests {
|
||||
|
||||
#[test]
|
||||
fn random_weights() {
|
||||
let mut rng = StdRng::from_entropy();
|
||||
let num_windows = 5;
|
||||
let weights = generate_random_weights(num_windows);
|
||||
let weights = generate_random_weights(num_windows, &mut rng);
|
||||
assert_eq!(num_windows, weights.len());
|
||||
for w in weights {
|
||||
assert!(0.0 < w && w < 1.0);
|
||||
@@ -109,8 +110,9 @@ mod tests {
|
||||
|
||||
#[test]
|
||||
fn random_weighted_dataset() {
|
||||
let mut rng = StdRng::from_entropy();
|
||||
let ds = build_hist_set();
|
||||
let rnd_weights_ds = generate_random_weighted_dataset(ds);
|
||||
let rnd_weights_ds = generate_random_weighted_dataset(ds, &mut rng);
|
||||
println!("{:?}", rnd_weights_ds.weights);
|
||||
for w in rnd_weights_ds.weights {
|
||||
assert!(w > 0.0);
|
||||
|
||||
394
src/io.rs
394
src/io.rs
@@ -1,6 +1,8 @@
|
||||
use super::histogram::Dataset;
|
||||
use super::histogram::Histogram;
|
||||
use super::Config;
|
||||
use super::correlation_analysis::{statistical_ineff, autocorrelation_time};
|
||||
use std::fs::OpenOptions;
|
||||
use std::fs::File;
|
||||
use std::io::prelude::*;
|
||||
use std::io::{BufReader,BufWriter};
|
||||
@@ -25,22 +27,26 @@ 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) -> Result<Dataset> {
|
||||
// given in the metadata file. This generates at least one Dataset,
|
||||
// or multiple Datasets if convdt is set in the config
|
||||
pub fn read_data(cfg: &Config) -> Result<Vec<Dataset>> {
|
||||
let mut bias_pos: Vec<f64> = Vec::new();
|
||||
let mut bias_fc: Vec<f64> = Vec::new();
|
||||
let mut histograms: Vec<Histogram> = Vec::new();
|
||||
let mut histograms: Vec<Vec<Histogram>> = Vec::new();
|
||||
let mut timeseries_lengths: Vec<usize> = Vec::new();
|
||||
let mut paths = Vec::new();
|
||||
|
||||
let kT = cfg.temperature * k_B;
|
||||
let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
|
||||
(cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
|
||||
}).collect();
|
||||
let num_bins = cfg.num_bins.iter().product();
|
||||
let num_bins: usize = cfg.num_bins.iter().product();
|
||||
let dimens_length = cfg.num_bins.clone();
|
||||
|
||||
let f = File::open(&cfg.metadata_file).chain_err(|| "Failed to open metadata file")?;
|
||||
let buf = BufReader::new(&f);
|
||||
|
||||
|
||||
// read each metadata file line and parse it
|
||||
for (line_num,l) in buf.lines().enumerate() {
|
||||
let line = l.chain_err(|| "Failed to read line")?;
|
||||
@@ -55,17 +61,6 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
|
||||
bail!(format!("Wrong number of columns in line {} of metadata file. Empty Line?", line_num+1));
|
||||
}
|
||||
|
||||
// parse histogram data
|
||||
let path = get_relative_path(&cfg.metadata_file, split[0]);
|
||||
let h = read_window_file(&path, cfg)
|
||||
.chain_err(|| format!("Failed to parse process data file {}", &path))?;
|
||||
if h.num_points == 0 {
|
||||
bail!(format!("No data points in histogram boundaries: {}", &path))
|
||||
}
|
||||
histograms.push(h);
|
||||
vprintln(format!("{}, {} data points added.", &path,
|
||||
histograms.last().unwrap().num_points), cfg.verbose);
|
||||
|
||||
// parse bias force constants and positions
|
||||
for val in split.iter().skip(1).take(cfg.dimens) {
|
||||
let pos = val.parse()
|
||||
@@ -77,15 +72,156 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
|
||||
.chain_err(|| format!("Failed to read bias fc in line {} of metadata file", line_num+1))?;
|
||||
bias_fc.push(fc);
|
||||
}
|
||||
|
||||
// parse histogram data
|
||||
let path = get_relative_path(&cfg.metadata_file, split[0]);
|
||||
paths.push(path.clone());
|
||||
let (timeseries, timeseries_initial_lengths) = read_window_file(&path, cfg)
|
||||
.chain_err(|| format!("Failed to read time series from {}", &path))?;
|
||||
timeseries_lengths.push(timeseries_initial_lengths);
|
||||
|
||||
|
||||
// for each timeseries, histograms are build for slices according to
|
||||
// start..convdt, start..2*convdt, ...
|
||||
histograms.push(Vec::new());
|
||||
let h_idx = histograms.len()-1;
|
||||
let convdt_stops = get_convdt_boundaries(×eries[0], &cfg);
|
||||
for (idx, interval) in convdt_stops.iter().enumerate() {
|
||||
// build histogram for slice start.._stop
|
||||
let (start, stop) = interval;
|
||||
let timeseries_mask: Vec<bool> = (0..timeseries[0].len()).map(|i| {
|
||||
is_in_time_boundaries(timeseries[0][i], *start, *stop)
|
||||
}).collect();
|
||||
let hist = build_histogram_from_timeseries(×eries, ×eries_mask, cfg);
|
||||
histograms[h_idx].push(hist);
|
||||
|
||||
if (cfg.convdt == 0.00) || idx+1 == convdt_stops.len() {
|
||||
vprintln(format!("{}, {} data points added.", &path,
|
||||
histograms[h_idx].last().unwrap().num_points), cfg.verbose);
|
||||
break
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if !histograms.is_empty() {
|
||||
Ok(Dataset::new(num_bins, dimens_length, bin_width, cfg.hist_min.clone(), cfg.hist_max.clone(), bias_pos, bias_fc, kT, histograms, cfg.cyclic))
|
||||
|
||||
// Datasets are created from histograms.
|
||||
// Empty histograms result in an error when its the final dataset,
|
||||
// and a warning otherwise.
|
||||
let num_datasets: usize = histograms.iter().map(|h| h.len()).max().unwrap();
|
||||
let dataset_boundaries: Vec<(f64, f64)> = (0..num_datasets).map(|idx| {
|
||||
(cfg.start, cfg.start+(idx as f64 + 1.0)*cfg.convdt) }
|
||||
).collect();
|
||||
vprintln(format!("Generating {} datasets from histograms.", num_datasets), cfg.verbose);
|
||||
let datasets: Vec<Dataset> = (0..num_datasets).map(|idx| {
|
||||
let mut dataset_histograms: Vec<Histogram> = Vec::with_capacity(histograms.len());
|
||||
for (hs, path) in histograms.iter().zip(&paths) {
|
||||
if hs.len() > idx {
|
||||
dataset_histograms.push(hs[idx].clone())
|
||||
} else {
|
||||
let warning = format!("No data points for interval {}-{} in histogram boundaries: {}.",
|
||||
dataset_boundaries[idx].0, dataset_boundaries[idx].1 ,&path);
|
||||
if !cfg.ignore_empty && idx+1 == num_datasets {
|
||||
bail!(warning + " This is the final dataset.");
|
||||
} else {
|
||||
eprintln!("{}", warning);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Ok(Dataset::new(num_bins, dimens_length.clone(), bin_width.clone(),
|
||||
cfg.hist_min.clone(), cfg.hist_max.clone(), bias_pos.clone(),
|
||||
bias_fc.clone(), kT, dataset_histograms, cfg.cyclic))
|
||||
}).collect::<Result<Vec<Dataset>>>().chain_err(|| "Failed to create datasets.")?;
|
||||
|
||||
if datasets.is_empty() {
|
||||
bail!("No datasets created.")
|
||||
} else if datasets[0].histograms.is_empty() {
|
||||
bail!("Dataset has no associated data points.")
|
||||
} else {
|
||||
bail!("Histogram has no datapoints.")
|
||||
if datasets.len() > 1 {
|
||||
println!("Datasets:");
|
||||
println!("Dataset\t\tTime interval\t\tWindows\t\tN_total");
|
||||
for (idx, dataset) in datasets.iter().enumerate() {
|
||||
let n: u32 = dataset.histograms.iter().map(|h| h.num_points).sum();
|
||||
let mut stop = cfg.start+cfg.convdt*(idx+1) as f64;
|
||||
if stop > cfg.end {
|
||||
stop = cfg.end;
|
||||
}
|
||||
println!("{:?}\t\t{:?}-{:?}\t\t{:?}\t\t{:?}", idx+1, cfg.start, stop, dataset.histograms.len(), n);
|
||||
}
|
||||
}
|
||||
|
||||
let histograms = &datasets.last().unwrap().histograms;
|
||||
if cfg.uncorr {
|
||||
println!("Timeseries Correlation:");
|
||||
println!("Window\t\tN\t\tN_uncorr\tN/N_uncorr");
|
||||
for (idx, (n, h)) in timeseries_lengths.iter().zip(histograms.iter()).enumerate() {
|
||||
println!("{:?}\t\t{:?}\t\t{:?}\t\t{:.2}",
|
||||
idx+1, n, h.num_points, h.num_points as f64 / *n as f64);
|
||||
}
|
||||
let total_n = timeseries_lengths.iter().sum::<usize>() as f64;
|
||||
let total_h = histograms.iter().map(|h| h.num_points).sum::<u32>() as f64;
|
||||
println!("\t\t\t\t\tTotal:\t{:.2}", total_h/total_n);
|
||||
}
|
||||
|
||||
Ok(datasets)
|
||||
}
|
||||
}
|
||||
|
||||
// builds a time boundaries for datasets from convdt, timeseries start and end
|
||||
fn get_convdt_boundaries(timeseries: &[f64], cfg: &Config) -> Vec<(f64, f64)> {
|
||||
let mut last_timestep = *timeseries.last().unwrap();
|
||||
if last_timestep > cfg.end {
|
||||
last_timestep = cfg.end;
|
||||
}
|
||||
let mut first_timestep = *timeseries.first().unwrap();
|
||||
if first_timestep < cfg.start {
|
||||
first_timestep = cfg.start;
|
||||
}
|
||||
if cfg.convdt == 0.0 {
|
||||
vec![(0.0, last_timestep)]
|
||||
} else {
|
||||
let intervals: usize = ((last_timestep - first_timestep) / cfg.convdt).ceil() as usize;
|
||||
(1..intervals+1).map(|i| {
|
||||
i as f64 * cfg.convdt + first_timestep
|
||||
}).map(|end| { (first_timestep, end) }).collect()
|
||||
}
|
||||
}
|
||||
|
||||
// build a histogram from a timeseries
|
||||
// mask is used to filter the timeseries for selected frames
|
||||
fn build_histogram_from_timeseries(timeseries: &[Vec<f64>], mask: &[bool],
|
||||
cfg: &Config) -> Histogram {
|
||||
|
||||
// total number of bins is the product of all dimensions length
|
||||
let total_bins = cfg.num_bins.iter().product();
|
||||
|
||||
// bin width for each dimension: (max-min)/bins
|
||||
let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
|
||||
(cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
|
||||
}).collect();
|
||||
|
||||
// build histogram for slice start..convdt_stop
|
||||
let mut hist = vec![0.0; total_bins];
|
||||
for i in (0..timeseries[0].len()).filter(|i| mask[*i]) {
|
||||
let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens+1];
|
||||
for j in 0..values.len() {
|
||||
values[j] = timeseries[j][i];
|
||||
}
|
||||
|
||||
if is_in_hist_boundaries(&values[1..], cfg) {
|
||||
let bin_indeces: Vec<usize> = (0..cfg.dimens).map(|dimen: usize| {
|
||||
let val = values[dimen+1];
|
||||
((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize
|
||||
}).collect();
|
||||
let index = flat_index(&bin_indeces, &cfg.num_bins);
|
||||
hist[index] += 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
let num_points: f64 = hist.iter().sum();
|
||||
Histogram::new(num_points as u32, hist)
|
||||
}
|
||||
|
||||
// transforms a multidimensional index into a one dimensional index
|
||||
// indeces: multidimensional indeces
|
||||
// lengths: length of the matrix in each dimension
|
||||
@@ -108,33 +244,57 @@ fn is_in_hist_boundaries(values: &[f64], cfg: &Config) -> bool {
|
||||
}
|
||||
|
||||
// returns true given time in inside the time boundaries defined by cfg
|
||||
fn is_in_time_boundaries(time: f64, cfg: &Config) -> bool {
|
||||
if cfg.start <= time && time <= cfg.end {
|
||||
fn is_in_time_boundaries(time: f64, start: f64, end: f64) -> bool {
|
||||
if start <= time && time <= end {
|
||||
return true
|
||||
}
|
||||
false
|
||||
}
|
||||
|
||||
// parse a time series file into a histogram
|
||||
fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
|
||||
let f = File::open(window_file)
|
||||
// parse a time series file
|
||||
fn read_window_file(window_file: &str, cfg: &Config) -> Result<(Vec<Vec<f64>>, usize)> {
|
||||
let mut timeseries: Vec<Vec<f64>> = read_timeseries(window_file, cfg)?;
|
||||
|
||||
// filter the timeseries based on start/end parameters
|
||||
let time_series_mask: Vec<bool> = timeseries[0].iter()
|
||||
.map(|t| is_in_time_boundaries(*t, cfg.start, cfg.end)).collect();
|
||||
timeseries = timeseries.into_iter().map(|ts| {
|
||||
ts.into_iter().zip(time_series_mask.iter()).filter_map(|(val, mask)| {
|
||||
if *mask {
|
||||
Some(val)
|
||||
} else {
|
||||
None
|
||||
}
|
||||
}).collect()
|
||||
}).collect::<Vec<Vec<f64>>>();
|
||||
|
||||
let timeseries_inital_length = timeseries[0].len();
|
||||
if cfg.uncorr {
|
||||
timeseries = uncorrelate(timeseries, cfg);
|
||||
}
|
||||
|
||||
if timeseries[0].is_empty() {
|
||||
bail!("Time series is empty")
|
||||
}
|
||||
|
||||
Ok((timeseries, timeseries_inital_length))
|
||||
}
|
||||
|
||||
// Read a multidimensional timeseries
|
||||
// The resulting vector contains one vector per dimension
|
||||
fn read_timeseries(window_file: &str, cfg: &Config) -> Result<Vec<Vec<f64>>> {
|
||||
let f = File::open(window_file)
|
||||
.chain_err(|| format!("Failed to open sample data file {}.", window_file))?;
|
||||
let mut buf = BufReader::new(&f);
|
||||
|
||||
// total number of bins is the product of all dimensions length
|
||||
let total_bins = cfg.num_bins.iter().product();
|
||||
let mut hist = vec![0.0; total_bins];
|
||||
|
||||
// bin width for each dimension: (max-min)/bins
|
||||
let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
|
||||
(cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
|
||||
}).collect();
|
||||
let mut timeseries = vec![Vec::new(); cfg.dimens+1];
|
||||
|
||||
// read and parse each timeseries line
|
||||
let mut line = String::new();
|
||||
let mut linecount = 0;
|
||||
while buf.read_line(&mut line).chain_err(|| "Failed to read line")? > 0 {
|
||||
linecount += 1;
|
||||
|
||||
// skip comments and empty lines
|
||||
if line.starts_with('#') || line.starts_with('@') || line.is_empty() {
|
||||
line.clear();
|
||||
@@ -142,43 +302,76 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
|
||||
}
|
||||
|
||||
{
|
||||
|
||||
let split: Vec<&str> = line.split_whitespace().collect();
|
||||
if split.len() < cfg.dimens+1 {
|
||||
bail!(format!("Wrong number of columns in line {} of window file {}. Empty Line?.", linecount, window_file));
|
||||
}
|
||||
|
||||
let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens+1];
|
||||
for i in 0..values.len() {
|
||||
values[i] = split[i].parse::<f64>()
|
||||
.chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))?;
|
||||
}
|
||||
for i in 0..cfg.dimens+1 {
|
||||
timeseries[i].push(split[i].parse::<f64>()
|
||||
.chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))?
|
||||
|
||||
if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) {
|
||||
let bin_indeces: Vec<usize> = (0..cfg.dimens).map(|dimen: usize| {
|
||||
let val = values[dimen+1];
|
||||
((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize
|
||||
}).collect();
|
||||
let index = flat_index(&bin_indeces, &cfg.num_bins);
|
||||
hist[index] += 1.0;
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
line.clear();
|
||||
}
|
||||
Ok(timeseries)
|
||||
}
|
||||
|
||||
let num_points: f64 = hist.iter().sum();
|
||||
Ok(Histogram::new(num_points as u32, hist))
|
||||
|
||||
// calculates the inefficiency for every collective variable
|
||||
// filters the timeseries based on the highest inefficiency
|
||||
fn uncorrelate(timeseries: Vec<Vec<f64>>, cfg: &Config) -> Vec<Vec<f64>> {
|
||||
// calculate inefficiencies and find the highest one
|
||||
let gs: Vec<f64> = timeseries[1..].iter().map(|ts| statistical_ineff(ts)).collect();
|
||||
let mut max_g = 1.0;
|
||||
for g in gs {
|
||||
if g > max_g {
|
||||
max_g = g;
|
||||
}
|
||||
}
|
||||
|
||||
// round g up
|
||||
let mut trunc_g = max_g.trunc() as usize;
|
||||
if (trunc_g as f64 - max_g).abs() > 0.000_000_000_1 {
|
||||
trunc_g += 1;
|
||||
}
|
||||
|
||||
// filter correlated samples from timeseries
|
||||
let prev_len = timeseries[0].len();
|
||||
let timeseries = timeseries.into_iter().map(|ts| {
|
||||
ts.into_iter().step_by(trunc_g).collect::<Vec<f64>>()
|
||||
}).collect::<Vec<Vec<f64>>>();
|
||||
|
||||
let new_len = timeseries[0].len();
|
||||
if cfg.verbose {
|
||||
let tau = autocorrelation_time(max_g)* (timeseries[0][1]-timeseries[0][0]);
|
||||
vprintln(format!("{:?}/{:?} samples are uncorrelated. {:?} samples removed from timeseries (tau={:.5})", new_len, prev_len, prev_len-new_len, tau), true);
|
||||
}
|
||||
timeseries
|
||||
}
|
||||
|
||||
// Write WHAM calculation results to out_file.
|
||||
pub fn write_results(out_file: &str, ds: &Dataset, free: &[f64],
|
||||
free_std: &[f64], prob: &[f64], prob_std: &[f64]) -> Result<()> {
|
||||
let output = File::create(out_file)
|
||||
pub fn write_results(out_file: &str, append: bool, ds: &Dataset, free: &[f64],
|
||||
free_std: &[f64], prob: &[f64], prob_std: &[f64], index: Option<usize>) -> Result<()> {
|
||||
|
||||
if !append && Path::new(out_file).exists() {
|
||||
std::fs::remove_file(out_file).chain_err(|| "Failed to delete file.")?;
|
||||
}
|
||||
let output = OpenOptions::new().write(true)
|
||||
.append(true)
|
||||
.create(true)
|
||||
.open(out_file)
|
||||
.chain_err(|| format!("Failed to create file with path {}", out_file))?;
|
||||
let mut buf = BufWriter::new(output);
|
||||
|
||||
let header: String = (0..ds.dimens_lengths.len()).map(|d| format!("coord{}", d+1))
|
||||
.collect::<Vec<String>>().join(" ");
|
||||
if let Some(index) = index {
|
||||
writeln!(buf, "#Dataset {}", index).unwrap();
|
||||
}
|
||||
writeln!(buf, "#{} Free Energy +/- Probability +/-", header).unwrap();
|
||||
|
||||
for bin in 0..free.len() {
|
||||
@@ -214,7 +407,10 @@ mod tests {
|
||||
bootstrap: 0,
|
||||
bootstrap_seed: 1234,
|
||||
start: 0.0,
|
||||
end: 1e+20
|
||||
end: 1e+20,
|
||||
uncorr: false,
|
||||
convdt: 0.0,
|
||||
ignore_empty: false,
|
||||
}
|
||||
}
|
||||
|
||||
@@ -222,8 +418,11 @@ mod tests {
|
||||
fn read_window_file() {
|
||||
let f = "example/1d_cyclic/COLVAR+0.0.xvg";
|
||||
let cfg = cfg();
|
||||
let h = super::read_window_file(&f, &cfg).unwrap();
|
||||
let (timeseries, timeseries_inital_length) = super::read_window_file(&f, &cfg).unwrap();
|
||||
let mask = vec![true; timeseries[0].len()];
|
||||
let h = build_histogram_from_timeseries(×eries, &mask, &cfg);
|
||||
println!("{:?}", h);
|
||||
assert_eq!(5000, timeseries_inital_length);
|
||||
assert_eq!(5000, h.num_points);
|
||||
assert_approx_eq!(0.0, h.bins[2]);
|
||||
assert_approx_eq!(11.0, h.bins[3]);
|
||||
@@ -233,11 +432,31 @@ mod tests {
|
||||
assert_approx_eq!(0.0, h.bins[7]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn read_timeseries() {
|
||||
let f = "example/1d_cyclic/COLVAR+0.0.xvg";
|
||||
let cfg = cfg();
|
||||
let ts = super::read_timeseries(&f, &cfg).unwrap();
|
||||
let expected = [
|
||||
-0.153_145,
|
||||
-0.377_860,
|
||||
0.010_992,
|
||||
0.123_074,
|
||||
0.108_291,
|
||||
0.261_607,
|
||||
];
|
||||
assert!(ts.len() == 2);
|
||||
assert!(ts[0].len() == 5000);
|
||||
println!("{:?}", ts);
|
||||
for (actual, expected) in ts[1].iter().zip(expected.iter()) {
|
||||
assert!((actual-expected).abs() < 0.001, format!("{:?} != {:?}", actual, expected));
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn read_data() {
|
||||
let cfg = cfg();
|
||||
let ds = super::read_data(&cfg).unwrap();
|
||||
let ds = &super::read_data(&cfg).unwrap()[0];
|
||||
println!("{:?}", ds);
|
||||
assert_eq!(25, ds.num_windows);
|
||||
assert_eq!(cfg.num_bins.len(), ds.dimens_lengths.len());
|
||||
@@ -256,4 +475,73 @@ mod tests {
|
||||
let relative3 = super::get_relative_path(&path1, &path3);
|
||||
assert_eq!("path/to/subfolder/another_file.dat" ,relative3);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn is_in_time_boundaries() {
|
||||
let start = 10.0;
|
||||
let end = 20.0;
|
||||
assert!(super::is_in_time_boundaries(15.0, start, end));
|
||||
assert!(super::is_in_time_boundaries(10.0, start, end));
|
||||
assert!(super::is_in_time_boundaries(20.0, start, end));
|
||||
assert!(!super::is_in_time_boundaries(9.9999999, start, end));
|
||||
assert!(!super::is_in_time_boundaries(20.000001, start, end));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn get_convdt_boundaries() {
|
||||
let mut cfg = cfg();
|
||||
|
||||
let timeseries: Vec<f64> = (0..31).map(|i| i as f64).collect();
|
||||
println!("{:?}", timeseries);
|
||||
|
||||
cfg.start = 10.0;
|
||||
cfg.end = 20.0;
|
||||
cfg.convdt = 10.0;
|
||||
let test = super::get_convdt_boundaries(×eries, &cfg);
|
||||
println!("{:?}", test);
|
||||
assert!(test.len() == 1);
|
||||
assert_approx_eq!(test[0].0, 10.0);
|
||||
assert_approx_eq!(test[0].1, 20.0);
|
||||
|
||||
cfg.start = 10.0;
|
||||
cfg.end = 20.0;
|
||||
cfg.convdt = 5.0;
|
||||
let test = super::get_convdt_boundaries(×eries, &cfg);
|
||||
println!("{:?}", test);
|
||||
assert!(test.len() == 2);
|
||||
assert_approx_eq!(test[0].0, 10.0);
|
||||
assert_approx_eq!(test[0].1, 15.0);
|
||||
assert_approx_eq!(test[1].0, 10.0);
|
||||
assert_approx_eq!(test[1].1, 20.0);
|
||||
|
||||
let timeseries: Vec<f64> = (10..21).map(|i| i as f64).collect();
|
||||
println!("{:?}", timeseries);
|
||||
|
||||
cfg.start = 10.0;
|
||||
cfg.end = 20.0;
|
||||
cfg.convdt = 10.0;
|
||||
let test = super::get_convdt_boundaries(×eries, &cfg);
|
||||
println!("{:?}", test);
|
||||
assert!(test.len() == 1);
|
||||
assert_approx_eq!(test[0].0, 10.0);
|
||||
assert_approx_eq!(test[0].1, 20.0);
|
||||
|
||||
cfg.start = 5.0;
|
||||
cfg.end = 20.0;
|
||||
cfg.convdt = 10.0;
|
||||
let test = super::get_convdt_boundaries(×eries, &cfg);
|
||||
println!("{:?}", test);
|
||||
assert!(test.len() == 1);
|
||||
assert_approx_eq!(test[0].0, 10.0);
|
||||
assert_approx_eq!(test[0].1, 20.0);
|
||||
|
||||
cfg.start = 5.0;
|
||||
cfg.end = 30.0;
|
||||
cfg.convdt = 10.0;
|
||||
let test = super::get_convdt_boundaries(×eries, &cfg);
|
||||
println!("{:?}", test);
|
||||
assert!(test.len() == 1);
|
||||
assert_approx_eq!(test[0].0, 10.0);
|
||||
assert_approx_eq!(test[0].1, 20.0);
|
||||
}
|
||||
}
|
||||
59
src/lib.rs
59
src/lib.rs
@@ -3,7 +3,6 @@
|
||||
#[macro_use]
|
||||
extern crate error_chain;
|
||||
extern crate rand;
|
||||
extern crate rgsl;
|
||||
extern crate rayon;
|
||||
#[cfg(test)]
|
||||
#[macro_use]
|
||||
@@ -13,6 +12,8 @@ extern crate assert_approx_eq;
|
||||
pub mod io;
|
||||
pub mod histogram;
|
||||
pub mod error_analysis;
|
||||
pub mod correlation_analysis;
|
||||
pub mod statistics;
|
||||
|
||||
use histogram::Dataset;
|
||||
use std::f64;
|
||||
@@ -45,16 +46,21 @@ pub struct Config {
|
||||
pub bootstrap_seed: u64,
|
||||
pub start: f64,
|
||||
pub end: f64,
|
||||
pub uncorr: bool,
|
||||
pub convdt: f64,
|
||||
pub ignore_empty: bool
|
||||
}
|
||||
|
||||
impl fmt::Display for Config {
|
||||
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
|
||||
write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?},
|
||||
verbose={}, tolerance={}, iterations={}, temperature={},
|
||||
cyclic={:?}, bootstrap={:?}, seed={:?}",
|
||||
cyclic={:?}, uncorr={:?}, bootstrap={:?}, seed={:?},
|
||||
uncorr={:?}, start={:?}, end={:?}, convdt={:?}, ignore_empty={:?}",
|
||||
self.metadata_file, self.hist_min, self.hist_max, self.num_bins,
|
||||
self.verbose, self.tolerance, self.max_iterations, self.temperature,
|
||||
self.cyclic, self.bootstrap, self.bootstrap_seed)
|
||||
self.cyclic, self.uncorr, self.bootstrap, self.bootstrap_seed,
|
||||
self.uncorr, self.start, self.end, self.convdt, self.ignore_empty)
|
||||
}
|
||||
}
|
||||
|
||||
@@ -181,27 +187,40 @@ pub fn run(cfg: &Config) -> Result<()>{
|
||||
println!("Supplied WHAM options: {}", &cfg);
|
||||
|
||||
println!("Reading input files.");
|
||||
let dataset = io::read_data(&cfg).chain_err(|| "Failed to create histogram.")?;
|
||||
println!("{}", &dataset);
|
||||
let datasets = io::read_data(&cfg).chain_err(|| "Failed to read data.")?;
|
||||
|
||||
let (P, F, F_prev) = perform_wham(&cfg, &dataset)?;
|
||||
println!("WHAM converged.");
|
||||
for (idx, dataset) in datasets.iter().enumerate() {
|
||||
if datasets.len() > 1 {
|
||||
println!("Dataset {}/{}: {}", idx+1, datasets.len(), &dataset);
|
||||
}
|
||||
else {
|
||||
println!("{}", &dataset);
|
||||
}
|
||||
let (P, F, F_prev) = perform_wham(&cfg, &dataset)?;
|
||||
println!("WHAM converged.");
|
||||
|
||||
let (P_std, free_energy_std) = if cfg.bootstrap > 0 {
|
||||
println!("Bootstrapping..");
|
||||
error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap)
|
||||
} else {
|
||||
(vec![0.0; P.len()], vec![0.0; P.len()])
|
||||
};
|
||||
let (P_std, free_energy_std) = if cfg.bootstrap > 0 {
|
||||
println!("Bootstrapping..");
|
||||
error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap)
|
||||
} else {
|
||||
(vec![0.0; P.len()], vec![0.0; P.len()])
|
||||
};
|
||||
|
||||
// calculate free energy and dump state
|
||||
println!("Finished. Dumping final PMF");
|
||||
let free_energy = calc_free_energy(&dataset, &P);
|
||||
dump_state(&dataset, &F, &F_prev, &P, &P_std, &free_energy, &free_energy_std);
|
||||
|
||||
io::write_results(&cfg.output, &dataset, &free_energy, &free_energy_std, &P, &P_std)
|
||||
.chain_err(|| "Could not write results to output file")?;
|
||||
// calculate free energy and dump state
|
||||
println!("Finished. Dumping PMF");
|
||||
let free_energy = calc_free_energy(&dataset, &P);
|
||||
|
||||
dump_state(&dataset, &F, &F_prev, &P, &P_std, &free_energy, &free_energy_std);
|
||||
let append = idx > 0 && datasets.len() > 1;
|
||||
let index = if datasets.len() > 1 {
|
||||
Some(idx)
|
||||
} else {
|
||||
None
|
||||
};
|
||||
io::write_results(&cfg.output, append, &dataset, &free_energy, &free_energy_std, &P, &P_std, index)
|
||||
.chain_err(|| "Could not write results to output file")?;
|
||||
}
|
||||
|
||||
Ok(())
|
||||
}
|
||||
|
||||
|
||||
@@ -58,6 +58,8 @@ fn cli() -> Result<Config> {
|
||||
.chain_err(|| "Cannot parse start time.")?;
|
||||
let end: f64 = matches.value_of("end").unwrap_or("1e+20").parse()
|
||||
.chain_err(|| "Cannot parse end time.")?;
|
||||
|
||||
let uncorr: bool = matches.is_present("uncorr");
|
||||
|
||||
if num_bins.len() != hist_max.len() || num_bins.len() != hist_max.len() {
|
||||
eprintln!("Input dimensions do not match (min: {}, max: {}, bins: {})",
|
||||
@@ -66,10 +68,14 @@ fn cli() -> Result<Config> {
|
||||
}
|
||||
|
||||
let dimens = num_bins.len();
|
||||
let convdt: f64 = matches.value_of("convdt").unwrap_or("0").parse()
|
||||
.chain_err(|| "Cannot parse convdt.")?;
|
||||
|
||||
let ignore_empty: bool = matches.is_present("ignore_empty");
|
||||
|
||||
Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens,
|
||||
verbose, tolerance, max_iterations, temperature, cyclic, output,
|
||||
bootstrap, bootstrap_seed, start, end})
|
||||
bootstrap, bootstrap_seed, start, end, uncorr, convdt, ignore_empty})
|
||||
}
|
||||
|
||||
fn main() {
|
||||
|
||||
59
src/statistics.rs
Normal file
59
src/statistics.rs
Normal file
@@ -0,0 +1,59 @@
|
||||
pub fn mean(x: &[f64]) -> f64 {
|
||||
x.iter().sum::<f64>() / x.len() as f64
|
||||
}
|
||||
|
||||
pub fn autocov(x: &[f64]) -> f64 {
|
||||
let x_mean = mean(x);
|
||||
|
||||
x.iter().map(|xi| {
|
||||
(xi-x_mean).powi(2)
|
||||
}).sum::<f64>() / x.len() as f64
|
||||
}
|
||||
|
||||
pub fn sd(x: &[f64]) -> f64 {
|
||||
let x_mean = mean(x);
|
||||
let n = x.len() as f64;
|
||||
|
||||
let sum = x.iter().map(|xi| {
|
||||
(xi-x_mean).powi(2)
|
||||
}).sum::<f64>();
|
||||
(1.0/(n-1.0) * sum).sqrt()
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use assert_approx_eq::assert_approx_eq;
|
||||
|
||||
// a sine wave
|
||||
fn dataset() -> Vec<f64> {
|
||||
(0..100).map(|i| (i as f64 / 100.0 * std::f64::consts::PI).sin())
|
||||
.collect::<Vec<f64>>()
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mean() {
|
||||
let ds = dataset();
|
||||
let expected =0.6366;
|
||||
|
||||
let m = super::mean(&ds);
|
||||
assert_approx_eq!(m, expected, 0.0001);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn autocorr() {
|
||||
let ds = dataset();
|
||||
let expected = 0.094_782;
|
||||
|
||||
let m = super::autocov(&ds);
|
||||
assert_approx_eq!(m, expected, 0.000_001);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sd() {
|
||||
let ds = dataset();
|
||||
let expected = 0.309_418;
|
||||
|
||||
let m = super::sd(&ds);
|
||||
assert_approx_eq!(m, expected, 0.000_001);
|
||||
}
|
||||
}
|
||||
@@ -6,44 +6,154 @@ mod integration {
|
||||
use std::process::Command;
|
||||
use std::fs;
|
||||
use super::command::get_command;
|
||||
use std::fs::OpenOptions;
|
||||
use std::io::prelude::*;
|
||||
|
||||
#[test]
|
||||
fn wham_1d_cyclic() {
|
||||
let output_file = "/tmp/wham_test_1d_cyclic.out";
|
||||
get_command()
|
||||
.args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"])
|
||||
.args(&["--bt", "100", "--seed", "1234"])
|
||||
.args(&["--seed", "1234"])
|
||||
.args(&["-f", "example/1d_cyclic/metadata.dat"])
|
||||
.args(&["-o", "/tmp/wham_test_1d_cyclic.out"])
|
||||
.args(&["-o", output_file])
|
||||
.output()
|
||||
.expect("failed to execute process");
|
||||
|
||||
assert!(fs::metadata("/tmp/wham_test_1d_cyclic.out").is_ok());
|
||||
assert!(fs::metadata(output_file).is_ok());
|
||||
let output = Command::new("diff")
|
||||
.arg("/tmp/wham_test_1d_cyclic.out")
|
||||
.arg(output_file)
|
||||
.arg("example/1d_cyclic/wham.out")
|
||||
.output()
|
||||
.expect("failed to run diff");
|
||||
let output_len = String::from_utf8_lossy(&output.stdout).len();
|
||||
assert_eq!(output_len, 0);
|
||||
std::fs::remove_file(output_file).unwrap();
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn wham_1d_cyclic_uncorrelated() {
|
||||
let output_file = "/tmp/wham_test_1d_cyclic.out";
|
||||
get_command()
|
||||
.args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic", "--uncorr"])
|
||||
.args(&["--seed", "1234"])
|
||||
.args(&["-f", "example/1d_cyclic/metadata.dat"])
|
||||
.args(&["-o", output_file])
|
||||
.output()
|
||||
.expect("failed to execute process");
|
||||
|
||||
assert!(fs::metadata("/tmp/wham_test_1d_cyclic.out").is_ok());
|
||||
let output = Command::new("diff")
|
||||
.arg(output_file)
|
||||
.arg("example/1d_cyclic/wham_uncorrelated.out")
|
||||
.output()
|
||||
.expect("failed to run diff");
|
||||
let output_len = String::from_utf8_lossy(&output.stdout).len();
|
||||
assert_eq!(output_len, 0);
|
||||
std::fs::remove_file(output_file).unwrap();
|
||||
}
|
||||
|
||||
#[test]
|
||||
// Test if convdt runs have the same result as normal runs
|
||||
fn wham_convdt() {
|
||||
// run wham with convdt
|
||||
let output_file = "/tmp/wham_test_convdt.out";
|
||||
get_command()
|
||||
.args(&["--bins", "10", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"])
|
||||
.args(&["--seed", "1234", "--tolerance", "0.001"])
|
||||
.args(&["--start", "0", "--end", "10"])
|
||||
.args(&["--convdt", "1"])
|
||||
.args(&["-f", "example/1d_cyclic/metadata.dat"])
|
||||
.args(&["-o", output_file])
|
||||
.output()
|
||||
.expect("failed to execute process");
|
||||
assert!(fs::metadata(output_file).is_ok());
|
||||
|
||||
// run wham for individual sets
|
||||
for i in 1..11 {
|
||||
let output_file_single = format!("/tmp/wham_test_convdt_{}.out", i);
|
||||
get_command()
|
||||
.args(&["--bins", "10", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"])
|
||||
.args(&["--seed", "1234", "--tolerance", "0.001"])
|
||||
.args(&["--start", "0", "--end", &i.to_string()])
|
||||
.args(&["-f", "example/1d_cyclic/metadata.dat"])
|
||||
.args(&["-o", &output_file_single])
|
||||
.output()
|
||||
.expect("failed to execute process");
|
||||
assert!(fs::metadata(output_file_single).is_ok());
|
||||
}
|
||||
|
||||
// combine individual runs
|
||||
let output_combined = "/tmp/wham_test_convdt_combined.out";
|
||||
let mut file = OpenOptions::new()
|
||||
.create(true)
|
||||
.write(true)
|
||||
.open(output_combined)
|
||||
.unwrap();
|
||||
for i in 1..11 {
|
||||
let output_file_single = format!("/tmp/wham_test_convdt_{}.out", i);
|
||||
println!("{}", output_file_single);
|
||||
file.write_all(format!("#Dataset {}\n", i-1).as_bytes()).unwrap();
|
||||
file.write_all(fs::read_to_string(output_file_single.clone()).unwrap().as_bytes()).unwrap();
|
||||
std::fs::remove_file(output_file_single).unwrap();
|
||||
}
|
||||
|
||||
// compare combined runs with single run
|
||||
let output = Command::new("diff")
|
||||
.arg(output_file)
|
||||
.arg(output_combined)
|
||||
.output()
|
||||
.expect("failed to run diff");
|
||||
let output_len = String::from_utf8_lossy(&output.stdout).len();
|
||||
assert_eq!(output_len, 0);
|
||||
|
||||
std::fs::remove_file(output_combined).unwrap();
|
||||
std::fs::remove_file(output_file).unwrap();
|
||||
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[ignore]
|
||||
fn wham_1d_cyclic_bootstrap() {
|
||||
let output_file = "/tmp/wham_test_1d_cyclic_bt.out";
|
||||
get_command()
|
||||
.args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"])
|
||||
.args(&["--seed", "1234", "--bt", "100"])
|
||||
.args(&["-f", "example/1d_cyclic/metadata.dat"])
|
||||
.args(&["-o", output_file])
|
||||
.output()
|
||||
.expect("failed to execute process");
|
||||
|
||||
assert!(fs::metadata(output_file).is_ok());
|
||||
let output = Command::new("diff")
|
||||
.arg(output_file)
|
||||
.arg("example/1d_cyclic/wham_bt.out")
|
||||
.output()
|
||||
.expect("failed to run diff");
|
||||
let output_len = String::from_utf8_lossy(&output.stdout).len();
|
||||
assert_eq!(output_len, 0);
|
||||
std::fs::remove_file(output_file).unwrap();
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[ignore] // expensive
|
||||
fn wham_2d_cyclic() {
|
||||
let output_file = "/tmp/wham_test_2d_cyclic.out";
|
||||
get_command()
|
||||
.args(&["--bins", "100,100", "--max", "pi,pi", "--min", "-pi,-pi", "-T", "300", "--cyclic"])
|
||||
.args(&["-f", "example/2d_cyclic/metadata.dat"])
|
||||
.args(&["-o", "/tmp/wham_test_2d_cyclic.out"])
|
||||
.args(&["-o", output_file])
|
||||
.output()
|
||||
.expect("failed to execute process");
|
||||
|
||||
assert!(fs::metadata("/tmp/wham_test_2d_cyclic.out").is_ok());
|
||||
assert!(fs::metadata(output_file).is_ok());
|
||||
let output = Command::new("diff")
|
||||
.arg("/tmp/wham_test_2d_cyclic.out")
|
||||
.arg(output_file)
|
||||
.arg("example/2d_cyclic/wham.out")
|
||||
.output()
|
||||
.expect("failed to run diff");
|
||||
let output_len = String::from_utf8_lossy(&output.stdout).len();
|
||||
assert_eq!(output_len, 0);
|
||||
std::fs::remove_file(output_file).unwrap();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user