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

Compare commits

base: bauer:v1.0.0
Branches Tags
bauer:main bauer:dependabot/cargo/rayon-1.10.0
bauer:v1.1.3 bauer:v1.1.2 bauer:v1.1.0 bauer:v1.0.0 bauer:v0.9.9 bauer:v0.9.8 bauer:0.9.8 bauer:v0.9.6 bauer:0.9.6 bauer:v0.9.3 bauer:v0.9.2 bauer:0.9.1 bauer:0.9.0
..
compare: bauer:v1.1.3
Branches Tags
bauer:dependabot/cargo/rayon-1.10.0 bauer:main
bauer:v1.1.3 bauer:v1.1.2 bauer:v1.1.0 bauer:v1.0.0 bauer:v0.9.9 bauer:v0.9.8 bauer:0.9.8 bauer:v0.9.6 bauer:0.9.6 bauer:v0.9.3 bauer:v0.9.2 bauer:0.9.1 bauer:0.9.0

15 Commits
v1.0.0 ... v1.1.3

Author SHA1 Message Date
Daniel Bauer
b6f338058a faster tests 2021-07-21 23:25:10 +02:00
Daniel Bauer
e0d2c1375a test convdt 2021-07-21 23:10:09 +02:00
Daniel Bauer
8802979a8e delete test files 2021-07-21 22:43:32 +02:00
Daniel Bauer
2d32d795f5 fix unit tests 2021-07-21 22:43:25 +02:00
Daniel Bauer
35ca1b1317 version up 2021-07-21 16:47:32 +02:00
Daniel Bauer
d9a219a36a better verbose output during dataset generation and flag to ignore empty histograms 2021-07-21 16:45:58 +02:00
Daniel Bauer
2c565dee28 remove some debug statements 2021-07-19 08:32:57 +02:00
Daniel Bauer
c05becea22 version up 2021-07-17 14:44:09 +02:00
Daniel Bauer
c27dbb7ab5 test g=1 for correlation 2021-07-17 14:42:21 +02:00
Daniel Bauer
c9e18ebda9 Squashed commit of the following: 
commit eaebf0dcbb259decbc0d8f5bbffa62244303f6c7
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 14:12:30 2021 +0200

    error for empty timeseries

commit 1fe5383c5079ca6c0ad102d4e2cb32d5b1227e80
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 14:09:13 2021 +0200

    refractored convdt slices calculation

commit 0b1fb7fb6a72edc50b013b59623551b2ccab6913
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 13:27:59 2021 +0200

    fix histogram building without convdt

commit f9882ca4cece57451cd2971d0993a2d3621b0cdf
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 12:53:06 2021 +0200

    fix tests not compiling

commit e6550e20bde3824f8199432b08777be41fc79fa8
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 12:43:36 2021 +0200

    refractoring

commit 18c77a9b6694d491ef23cebb3918ccda8fcc44a5
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Fri Jul 16 08:16:34 2021 +0200

    run and output for multiple datasets

commit 069f318f72207c416387007435eeff9867494e7a
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Fri Jul 16 07:58:01 2021 +0200

    cleanup io.rs

commit 10efa428a0c6bf490c9f2d7a4c1df185de402a11
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Thu Jul 15 19:27:22 2021 +0200

    parse multiple datasets with convdt
 2021-07-17 14:14:25 +02:00
Daniel Bauer
0b5d0ebecd test time boundaries 2021-07-17 13:31:36 +02:00
Daniel Bauer
091b1f1382 TODO 2021-07-17 13:04:42 +02:00
Daniel Bauer
9849c00321 cargo ignore examples 2021-07-17 11:30:42 +02:00
daniel
034586d08e inlined autocorrelation calculation for better performance 2020-10-26 11:24:48 +01:00
Daniel Bauer
cc4448f1d8 Update README.md 2020-10-26 10:53:57 +01:00
11 changed files with 633 additions and 241 deletions
Expand all files Collapse all files

2
Cargo.lock generated
View File

@@ -362,7 +362,7 @@ checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519"
[[package]] [[package]]
name = "wham" name = "wham"
version = "1.0.0" version = "1.1.2"
dependencies = [ dependencies = [
"assert_approx_eq", "assert_approx_eq",
"clap", "clap",

5
Cargo.toml
View File

@@ -1,6 +1,6 @@
[package] [package]
name = "wham" name = "wham"
version = "1.0.0" version = "1.1.2"
authors = ["Daniel Bauer <bauer@cbs.tu-darmstadt.de>"] authors = ["Daniel Bauer <bauer@cbs.tu-darmstadt.de>"]
description = "An implementation of the weighted histogram analysis method" description = "An implementation of the weighted histogram analysis method"
license = "GPL-3.0" license = "GPL-3.0"
@@ -8,6 +8,9 @@ repository = "https://github.com/danijoo/WHAM"
readme = "README.md" readme = "README.md"
categories = ["science", "command-line-utilities", "algorithms"] categories = ["science", "command-line-utilities", "algorithms"]
keywords = ["math", "statistics", "histogram", "bioinformatics", "molecular-dynamics"] keywords = ["math", "statistics", "histogram", "bioinformatics", "molecular-dynamics"]
exclude = [
"example/*"
]
[dependencies] [dependencies]
clap = {version="2.32.0", features=['yaml']} clap = {version="2.32.0", features=['yaml']}

17
README.md
View File

@@ -20,12 +20,6 @@ Features
Installation Installation
--- ---
WHAM requires the GSL library to be installed:
```bash
# on debian/ubuntu:
sudo apt-get install libgsl0-dev
```
Installation from source via cargo: Installation from source via cargo:
```bash ```bash
# cargo installation # cargo installation
@@ -40,8 +34,8 @@ wham has a convenient command line interface. You can see all options with
```wham -h```: ```wham -h```:
``` ```
wham 1.0.0 wham 1.1.0
D. Bauer <bauer@bio.tu-darmstadt.de> 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 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. potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
@@ -62,7 +56,7 @@ The first column will be ignored and is followed by N reaction coordinates x.
Shipped under the GPLv3 license. Shipped under the GPLv3 license.
USAGE: USAGE:
wham [FLAGS] [OPTIONS] --bins <BINS> --max <HIST_MAX> --file <METADATA> --min <HIST_MIN> --temperature <temperature> wham [FLAGS] [OPTIONS] --bins <BINS> --max <HIST_MAX> --file <METADATA> --min <HIST_MIN> --temperature <temperature>
FLAGS: FLAGS:
-c, --cyclic For periodic reaction coordinates. If this is set, the first and last coordinate bin in each -c, --cyclic For periodic reaction coordinates. If this is set, the first and last coordinate bin in each
@@ -78,6 +72,10 @@ OPTIONS:
--bt <bootstrap> Number of bayesian bootstrapping runs for error analysis by assigning random --bt <bootstrap> Number of bayesian bootstrapping runs for error analysis by assigning random
weights (defaults to 0). weights (defaults to 0).
--seed <bootstrap_seed> Random seed for bootstrapping runs. --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) --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). -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". --max <HIST_MAX> Histogram maxima (comma separated). Also accepts "pi".
@@ -158,6 +156,7 @@ tempering simulations, JCTC 3(1):26-41*
TODO TODO
--- ---
- Option to output histograms
- Replica exchange - Replica exchange
License & Citing License & Citing

200
example/1d_cyclic/wham.out
View File

@@ -1,101 +1,101 @@
#coord1 Free Energy +/- Probability +/- #coord1 Free Energy +/- Probability +/-
-3.110177 7.158102 0.071438 0.003494 0.000068 -3.110177 7.158102 0.000000 0.003494 0.000000
-3.047345 5.365727 0.069449 0.007168 0.000135 -3.047345 5.365727 0.000000 0.007168 0.000000
-2.984513 3.873190 0.067495 0.013039 0.000232 -2.984513 3.873190 0.000000 0.013039 0.000000
-2.921681 2.953162 0.067589 0.018855 0.000343 -2.921681 2.953162 0.000000 0.018855 0.000000
-2.858849 1.949554 0.064640 0.028195 0.000480 -2.858849 1.949554 0.000000 0.028195 0.000000
-2.796017 1.391747 0.063570 0.035261 0.000584 -2.796017 1.391747 0.000000 0.035261 0.000000
-2.733186 1.128270 0.061710 0.039189 0.000620 -2.733186 1.128270 0.000000 0.039189 0.000000
-2.670354 0.839970 0.060445 0.043991 0.000667 -2.670354 0.839970 0.000000 0.043991 0.000000
-2.607522 0.624769 0.060762 0.047955 0.000739 -2.607522 0.624769 0.000000 0.047955 0.000000
-2.544690 0.663757 0.060786 0.047211 0.000731 -2.544690 0.663757 0.000000 0.047211 0.000000
-2.481858 1.051932 0.059774 0.040407 0.000617 -2.481858 1.051932 0.000000 0.040407 0.000000
-2.419026 1.463048 0.060422 0.034267 0.000527 -2.419026 1.463048 0.000000 0.034267 0.000000
-2.356194 1.990616 0.055242 0.027734 0.000360 -2.356194 1.990616 0.000000 0.027734 0.000000
-2.293363 2.190692 0.044458 0.025597 0.000210 -2.293363 2.190692 0.000000 0.025597 0.000000
-2.230531 2.553036 0.042144 0.022136 0.000159 -2.230531 2.553036 0.000000 0.022136 0.000000
-2.167699 2.572522 0.043313 0.021964 0.000165 -2.167699 2.572522 0.000000 0.021964 0.000000
-2.104867 2.472360 0.039799 0.022863 0.000157 -2.104867 2.472360 0.000000 0.022863 0.000000
-2.042035 2.517562 0.036792 0.022453 0.000203 -2.042035 2.517562 0.000000 0.022453 0.000000
-1.979203 2.469778 0.036324 0.022887 0.000237 -1.979203 2.469778 0.000000 0.022887 0.000000
-1.916372 2.223125 0.049496 0.025266 0.000565 -1.916372 2.223125 0.000000 0.025266 0.000000
-1.853540 2.080157 0.039205 0.026756 0.000453 -1.853540 2.080157 0.000000 0.026756 0.000000
-1.790708 1.793841 0.041488 0.030011 0.000558 -1.790708 1.793841 0.000000 0.030011 0.000000
-1.727876 1.458784 0.042362 0.034326 0.000565 -1.727876 1.458784 0.000000 0.034326 0.000000
-1.665044 0.914441 0.037561 0.042697 0.000682 -1.665044 0.914441 0.000000 0.042697 0.000000
-1.602212 0.303813 0.028701 0.054540 0.000761 -1.602212 0.303813 0.000000 0.054540 0.000000
-1.539380 0.268335 0.025765 0.055321 0.000771 -1.539380 0.268335 0.000000 0.055321 0.000000
-1.476549 0.000000 0.023371 0.061604 0.000866 -1.476549 0.000000 0.000000 0.061604 0.000000
-1.413717 0.537061 0.026370 0.049671 0.000715 -1.413717 0.537061 0.000000 0.049671 0.000000
-1.350885 1.439940 0.025620 0.034586 0.000547 -1.350885 1.439940 0.000000 0.034586 0.000000
-1.288053 2.391838 0.027653 0.023614 0.000393 -1.288053 2.391838 0.000000 0.023614 0.000000
-1.225221 3.779470 0.027093 0.013538 0.000229 -1.225221 3.779470 0.000000 0.013538 0.000000
-1.162389 5.685555 0.027762 0.006305 0.000111 -1.162389 5.685555 0.000000 0.006305 0.000000
-1.099557 7.661896 0.028996 0.002855 0.000051 -1.099557 7.661896 0.000000 0.002855 0.000000
-1.036726 9.946594 0.032118 0.001142 0.000021 -1.036726 9.946594 0.000000 0.001142 0.000000
-0.973894 12.359408 0.042801 0.000434 0.000009 -0.973894 12.359408 0.000000 0.000434 0.000000
-0.911062 14.954655 0.052443 0.000153 0.000004 -0.911062 14.954655 0.000000 0.000153 0.000000
-0.848230 17.742242 0.063261 0.000050 0.000002 -0.848230 17.742242 0.000000 0.000050 0.000000
-0.785398 20.555785 0.068133 0.000016 0.000001 -0.785398 20.555785 0.000000 0.000016 0.000000
-0.722566 22.811160 0.075518 0.000007 0.000000 -0.722566 22.811160 0.000000 0.000007 0.000000
-0.659734 25.178094 0.086226 0.000003 0.000000 -0.659734 25.178094 0.000000 0.000003 0.000000
-0.596903 26.442288 0.091783 0.000002 0.000000 -0.596903 26.442288 0.000000 0.000002 0.000000
-0.534071 27.897565 0.093175 0.000001 0.000000 -0.534071 27.897565 0.000000 0.000001 0.000000
-0.471239 29.062473 0.096027 0.000001 0.000000 -0.471239 29.062473 0.000000 0.000001 0.000000
-0.408407 30.384521 0.096855 0.000000 0.000000 -0.408407 30.384521 0.000000 0.000000 0.000000
-0.345575 31.638454 0.099956 0.000000 0.000000 -0.345575 31.638454 0.000000 0.000000 0.000000
-0.282743 32.817727 0.107689 0.000000 0.000000 -0.282743 32.817727 0.000000 0.000000 0.000000
-0.219911 33.770369 0.108946 0.000000 0.000000 -0.219911 33.770369 0.000000 0.000000 0.000000
-0.157080 34.505503 0.112902 0.000000 0.000000 -0.157080 34.505503 0.000000 0.000000 0.000000
-0.094248 35.431659 0.120856 0.000000 0.000000 -0.094248 35.431659 0.000000 0.000000 0.000000
-0.031416 35.615810 0.122046 0.000000 0.000000 -0.031416 35.615810 0.000000 0.000000 0.000000
0.031416 35.561946 0.118225 0.000000 0.000000 0.031416 35.561946 0.000000 0.000000 0.000000
0.094248 35.382089 0.108652 0.000000 0.000000 0.094248 35.382089 0.000000 0.000000 0.000000
0.157080 34.934827 0.118384 0.000000 0.000000 0.157080 34.934827 0.000000 0.000000 0.000000
0.219911 33.673460 0.119195 0.000000 0.000000 0.219911 33.673460 0.000000 0.000000 0.000000
0.282743 32.731563 0.121592 0.000000 0.000000 0.282743 32.731563 0.000000 0.000000 0.000000
0.345575 31.261855 0.129698 0.000000 0.000000 0.345575 31.261855 0.000000 0.000000 0.000000
0.408407 29.717377 0.139789 0.000000 0.000000 0.408407 29.717377 0.000000 0.000000 0.000000
0.471239 28.076620 0.137289 0.000001 0.000000 0.471239 28.076620 0.000000 0.000001 0.000000
0.534071 26.481097 0.135712 0.000002 0.000000 0.534071 26.481097 0.000000 0.000002 0.000000
0.596903 24.487358 0.135156 0.000003 0.000000 0.596903 24.487358 0.000000 0.000003 0.000000
0.659734 22.344251 0.131027 0.000008 0.000000 0.659734 22.344251 0.000000 0.000008 0.000000
0.722566 20.241543 0.133145 0.000018 0.000001 0.722566 20.241543 0.000000 0.000018 0.000000
0.785398 18.341869 0.131952 0.000039 0.000002 0.785398 18.341869 0.000000 0.000039 0.000000
0.848230 16.261582 0.134676 0.000091 0.000006 0.848230 16.261582 0.000000 0.000091 0.000000
0.911062 14.301801 0.134417 0.000199 0.000013 0.911062 14.301801 0.000000 0.000199 0.000000
0.973894 12.603788 0.131018 0.000394 0.000024 0.973894 12.603788 0.000000 0.000394 0.000000
1.036726 11.249601 0.131276 0.000678 0.000043 1.036726 11.249601 0.000000 0.000678 0.000000
1.099557 10.087886 0.132498 0.001079 0.000070 1.099557 10.087886 0.000000 0.001079 0.000000
1.162389 9.443303 0.132990 0.001398 0.000092 1.162389 9.443303 0.000000 0.001398 0.000000
1.225221 9.152799 0.132250 0.001570 0.000104 1.225221 9.152799 0.000000 0.001570 0.000000
1.288053 9.331937 0.133178 0.001462 0.000099 1.288053 9.331937 0.000000 0.001462 0.000000
1.350885 9.905546 0.133357 0.001161 0.000078 1.350885 9.905546 0.000000 0.001161 0.000000
1.413717 11.042050 0.133807 0.000736 0.000051 1.413717 11.042050 0.000000 0.000736 0.000000
1.476549 12.598167 0.132722 0.000395 0.000027 1.476549 12.598167 0.000000 0.000395 0.000000
1.539380 14.520167 0.131816 0.000183 0.000012 1.539380 14.520167 0.000000 0.000183 0.000000
1.602212 16.569783 0.131773 0.000080 0.000005 1.602212 16.569783 0.000000 0.000080 0.000000
1.665044 18.687390 0.132601 0.000034 0.000002 1.665044 18.687390 0.000000 0.000034 0.000000
1.727876 20.775408 0.133867 0.000015 0.000001 1.727876 20.775408 0.000000 0.000015 0.000000
1.790708 22.905200 0.129265 0.000006 0.000000 1.790708 22.905200 0.000000 0.000006 0.000000
1.853540 24.643852 0.128894 0.000003 0.000000 1.853540 24.643852 0.000000 0.000003 0.000000
1.916372 26.301740 0.130266 0.000002 0.000000 1.916372 26.301740 0.000000 0.000002 0.000000
1.979203 27.372071 0.128620 0.000001 0.000000 1.979203 27.372071 0.000000 0.000001 0.000000
2.042035 28.697726 0.133263 0.000001 0.000000 2.042035 28.697726 0.000000 0.000001 0.000000
2.104867 29.417901 0.133513 0.000000 0.000000 2.104867 29.417901 0.000000 0.000000 0.000000
2.167699 30.008351 0.130925 0.000000 0.000000 2.167699 30.008351 0.000000 0.000000 0.000000
2.230531 30.406016 0.124398 0.000000 0.000000 2.230531 30.406016 0.000000 0.000000 0.000000
2.293363 30.171275 0.122493 0.000000 0.000000 2.293363 30.171275 0.000000 0.000000 0.000000
2.356194 29.884646 0.130278 0.000000 0.000000 2.356194 29.884646 0.000000 0.000000 0.000000
2.419026 29.428153 0.130844 0.000000 0.000000 2.419026 29.428153 0.000000 0.000000 0.000000
2.481858 28.546982 0.148114 0.000001 0.000000 2.481858 28.546982 0.000000 0.000001 0.000000
2.544690 27.757520 0.133045 0.000001 0.000000 2.544690 27.757520 0.000000 0.000001 0.000000
2.607522 26.505787 0.134364 0.000001 0.000000 2.607522 26.505787 0.000000 0.000001 0.000000
2.670354 24.491866 0.115061 0.000003 0.000000 2.670354 24.491866 0.000000 0.000003 0.000000
2.733186 22.320664 0.110036 0.000008 0.000000 2.733186 22.320664 0.000000 0.000008 0.000000
2.796017 20.052723 0.107894 0.000020 0.000001 2.796017 20.052723 0.000000 0.000020 0.000000
2.858849 17.655650 0.105964 0.000052 0.000002 2.858849 17.655650 0.000000 0.000052 0.000000
2.921681 15.471590 0.107005 0.000125 0.000005 2.921681 15.471590 0.000000 0.000125 0.000000
2.984513 13.138167 0.099378 0.000318 0.000012 2.984513 13.138167 0.000000 0.000318 0.000000
3.047345 11.092386 0.087770 0.000722 0.000022 3.047345 11.092386 0.000000 0.000722 0.000000
3.110177 9.065722 0.077092 0.001626 0.000038 3.110177 9.065722 0.000000 0.001626 0.000000

101
example/1d_cyclic/wham_bt.out Normal file
View 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

12
src/cli.yml
View File

@@ -1,5 +1,5 @@
name: wham name: wham
version: "1.0.0" version: "1.1.2"
author: D. Bauer <bauer@cbs.tu-darmstadt.de> author: D. Bauer <bauer@cbs.tu-darmstadt.de>
about: | 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. 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.
@@ -107,3 +107,13 @@ args:
help: Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated samples (default is off). help: Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated samples (default is off).
takes_value: false takes_value: false
required: 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

48
src/correlation_analysis.rs
View File

@@ -8,11 +8,14 @@ use super::statistics;
// Method for the Analysis of Simulated and Parallel Tempering Simulations, JCTC" // Method for the Analysis of Simulated and Parallel Tempering Simulations, JCTC"
pub fn statistical_ineff(timeseries: &[f64]) -> f64 { pub fn statistical_ineff(timeseries: &[f64]) -> f64 {
let n = timeseries.len(); let n = timeseries.len();
let autocorr = autocorrelation(timeseries); 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; let mut g = 1.0;
for t in 1..(n-1) { for t in 1..(n-1) {
let c = autocorr[t-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 { if c <= 0.0 {
break; break;
} }
@@ -25,22 +28,6 @@ pub fn statistical_ineff(timeseries: &[f64]) -> f64 {
} }
} }
// calculates the autocorrelation of a simeseries
fn autocorrelation(timeseries: &[f64]) -> Vec<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 autocorr = Vec::new();
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: f64 = tmp.map(|x| x+x).sum::<f64>() / (2.0 * (n as f64-t as f64)*cov);
autocorr.push(c);
}
autocorr
}
// The autocorrelation time of a timeseries can be deduced from the // The autocorrelation time of a timeseries can be deduced from the
// `statistical_ineff` by (g-1)/2.0 // `statistical_ineff` by (g-1)/2.0
pub fn autocorrelation_time(g: f64) -> f64 { pub fn autocorrelation_time(g: f64) -> f64 {
@@ -65,29 +52,18 @@ mod tests {
} }
#[test]
fn autocorrelation() {
let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
let autocorr = super::autocorrelation(&timeseries);
let expected = [
0.691_900_865_597_914_3, 0.533_171_939_967_135_5,
0.204_726_209_564_635_89, -0.002_850_876_458_920_514,
-0.138_428_500_779_381_46, -0.265_292_355_297_323_2,
-0.314_271_982_722_353_85, -0.261_750_515_155_769_3,
-0.205_948_647_302_903_38, -0.133_100_190_918_118_12,
-0.188_756_890_119_342_6, -0.194_493_662_542_493_3,
-0.196_359_967_318_946_1, -0.113_915_878_338_380_03];
for (actual, expected) in autocorr.iter().zip(expected.iter()) {
assert!((actual-expected).abs() < 0.001);
}
}
#[test] #[test]
fn statistical_ineff() { fn statistical_ineff() {
let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg"); let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
let g = super::statistical_ineff(&timeseries); let g = super::statistical_ineff(&timeseries);
println!("{:?}", g); println!("{:?}", g);
assert!((g - 3.859).abs() < 0.001) 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(&timeseries);
assert_approx_eq!(g, 1.0);
} }
#[test] #[test]

310
src/io.rs
View File

@@ -2,6 +2,7 @@ use super::histogram::Dataset;
use super::histogram::Histogram; use super::histogram::Histogram;
use super::Config; use super::Config;
use super::correlation_analysis::{statistical_ineff, autocorrelation_time}; use super::correlation_analysis::{statistical_ineff, autocorrelation_time};
use std::fs::OpenOptions;
use std::fs::File; use std::fs::File;
use std::io::prelude::*; use std::io::prelude::*;
use std::io::{BufReader,BufWriter}; use std::io::{BufReader,BufWriter};
@@ -26,23 +27,26 @@ pub fn vprintln(s: String, verbose: bool) {
} }
// Read input data into a histogram set by iterating over input files // Read input data into a histogram set by iterating over input files
// given in the metadata file // given in the metadata file. This generates at least one Dataset,
pub fn read_data(cfg: &Config) -> Result<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_pos: Vec<f64> = Vec::new();
let mut bias_fc: Vec<f64> = Vec::new(); let mut bias_fc: Vec<f64> = Vec::new();
let mut histograms: Vec<Histogram> = Vec::new(); let mut histograms: Vec<Vec<Histogram>> = Vec::new();
let mut timeseries_lengths: Vec<usize> = Vec::new(); let mut timeseries_lengths: Vec<usize> = Vec::new();
let mut paths = Vec::new();
let kT = cfg.temperature * k_B; let kT = cfg.temperature * k_B;
let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| { let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
(cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64) (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
}).collect(); }).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 dimens_length = cfg.num_bins.clone();
let f = File::open(&cfg.metadata_file).chain_err(|| "Failed to open metadata file")?; let f = File::open(&cfg.metadata_file).chain_err(|| "Failed to open metadata file")?;
let buf = BufReader::new(&f); let buf = BufReader::new(&f);
// read each metadata file line and parse it // read each metadata file line and parse it
for (line_num,l) in buf.lines().enumerate() { for (line_num,l) in buf.lines().enumerate() {
let line = l.chain_err(|| "Failed to read line")?; let line = l.chain_err(|| "Failed to read line")?;
@@ -57,18 +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)); 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, timeseries_inital_length) = 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);
timeseries_lengths.push(timeseries_inital_length);
vprintln(format!("{}, {} data points added.", &path,
histograms.last().unwrap().num_points), cfg.verbose);
// parse bias force constants and positions // parse bias force constants and positions
for val in split.iter().skip(1).take(cfg.dimens) { for val in split.iter().skip(1).take(cfg.dimens) {
let pos = val.parse() let pos = val.parse()
@@ -80,12 +72,87 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
.chain_err(|| format!("Failed to read bias fc in line {} of metadata file", line_num+1))?; .chain_err(|| format!("Failed to read bias fc in line {} of metadata file", line_num+1))?;
bias_fc.push(fc); 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(&timeseries[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(&timeseries, &timeseries_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() { // 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 {
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 { if cfg.uncorr {
println!("Timeseries Correlation"); println!("Timeseries Correlation:");
println!();
println!("Window\t\tN\t\tN_uncorr\tN/N_uncorr"); println!("Window\t\tN\t\tN_uncorr\tN/N_uncorr");
for (idx, (n, h)) in timeseries_lengths.iter().zip(histograms.iter()).enumerate() { for (idx, (n, h)) in timeseries_lengths.iter().zip(histograms.iter()).enumerate() {
println!("{:?}\t\t{:?}\t\t{:?}\t\t{:.2}", println!("{:?}\t\t{:?}\t\t{:?}\t\t{:.2}",
@@ -94,15 +161,67 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
let total_n = timeseries_lengths.iter().sum::<usize>() 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; 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); println!("\t\t\t\t\tTotal:\t{:.2}", total_h/total_n);
} }
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)) Ok(datasets)
} else {
bail!("Histogram has no datapoints.")
} }
} }
// 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 // transforms a multidimensional index into a one dimensional index
// indeces: multidimensional indeces // indeces: multidimensional indeces
// lengths: length of the matrix in each dimension // lengths: length of the matrix in each dimension
@@ -125,50 +244,40 @@ fn is_in_hist_boundaries(values: &[f64], cfg: &Config) -> bool {
} }
// returns true given time in inside the time boundaries defined by cfg // returns true given time in inside the time boundaries defined by cfg
fn is_in_time_boundaries(time: f64, cfg: &Config) -> bool { fn is_in_time_boundaries(time: f64, start: f64, end: f64) -> bool {
if cfg.start <= time && time <= cfg.end { if start <= time && time <= end {
return true return true
} }
false false
} }
// parse a time series file
// parse a time series file into a histogram fn read_window_file(window_file: &str, cfg: &Config) -> Result<(Vec<Vec<f64>>, usize)> {
fn read_window_file(window_file: &str, cfg: &Config) -> Result<(Histogram, usize)> {
// 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<f64>> = read_timeseries(window_file, cfg)?; let mut timeseries: Vec<Vec<f64>> = read_timeseries(window_file, cfg)?;
let timeseries_inital_length = timeseries[0].len();
// 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 { if cfg.uncorr {
timeseries = uncorrelate(timeseries, cfg); timeseries = uncorrelate(timeseries, cfg);
} }
for i in 0..timeseries[0].len() { if timeseries[0].is_empty() {
let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens+1]; bail!("Time series is empty")
for j in 0..values.len() {
values[j] = timeseries[j][i];
}
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;
}
} }
let num_points: f64 = hist.iter().sum(); Ok((timeseries, timeseries_inital_length))
Ok((Histogram::new(num_points as u32, hist), timeseries_inital_length))
} }
// Read a multidimensional timeseries // Read a multidimensional timeseries
@@ -245,14 +354,24 @@ fn uncorrelate(timeseries: Vec<Vec<f64>>, cfg: &Config) -> Vec<Vec<f64>> {
} }
// Write WHAM calculation results to out_file. // Write WHAM calculation results to out_file.
pub fn write_results(out_file: &str, ds: &Dataset, free: &[f64], pub fn write_results(out_file: &str, append: bool, ds: &Dataset, free: &[f64],
free_std: &[f64], prob: &[f64], prob_std: &[f64]) -> Result<()> { free_std: &[f64], prob: &[f64], prob_std: &[f64], index: Option<usize>) -> Result<()> {
let output = File::create(out_file)
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))?; .chain_err(|| format!("Failed to create file with path {}", out_file))?;
let mut buf = BufWriter::new(output); let mut buf = BufWriter::new(output);
let header: String = (0..ds.dimens_lengths.len()).map(|d| format!("coord{}", d+1)) let header: String = (0..ds.dimens_lengths.len()).map(|d| format!("coord{}", d+1))
.collect::<Vec<String>>().join(" "); .collect::<Vec<String>>().join(" ");
if let Some(index) = index {
writeln!(buf, "#Dataset {}", index).unwrap();
}
writeln!(buf, "#{} Free Energy +/- Probability +/-", header).unwrap(); writeln!(buf, "#{} Free Energy +/- Probability +/-", header).unwrap();
for bin in 0..free.len() { for bin in 0..free.len() {
@@ -290,6 +409,8 @@ mod tests {
start: 0.0, start: 0.0,
end: 1e+20, end: 1e+20,
uncorr: false, uncorr: false,
convdt: 0.0,
ignore_empty: false,
} }
} }
@@ -297,7 +418,9 @@ mod tests {
fn read_window_file() { fn read_window_file() {
let f = "example/1d_cyclic/COLVAR+0.0.xvg"; let f = "example/1d_cyclic/COLVAR+0.0.xvg";
let cfg = cfg(); let cfg = cfg();
let (h, timeseries_inital_length) = 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(&timeseries, &mask, &cfg);
println!("{:?}", h); println!("{:?}", h);
assert_eq!(5000, timeseries_inital_length); assert_eq!(5000, timeseries_inital_length);
assert_eq!(5000, h.num_points); assert_eq!(5000, h.num_points);
@@ -333,7 +456,7 @@ mod tests {
#[test] #[test]
fn read_data() { fn read_data() {
let cfg = cfg(); let cfg = cfg();
let ds = super::read_data(&cfg).unwrap(); let ds = &super::read_data(&cfg).unwrap()[0];
println!("{:?}", ds); println!("{:?}", ds);
assert_eq!(25, ds.num_windows); assert_eq!(25, ds.num_windows);
assert_eq!(cfg.num_bins.len(), ds.dimens_lengths.len()); assert_eq!(cfg.num_bins.len(), ds.dimens_lengths.len());
@@ -352,4 +475,73 @@ mod tests {
let relative3 = super::get_relative_path(&path1, &path3); let relative3 = super::get_relative_path(&path1, &path3);
assert_eq!("path/to/subfolder/another_file.dat" ,relative3); 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(&timeseries, &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(&timeseries, &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(&timeseries, &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(&timeseries, &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(&timeseries, &cfg);
println!("{:?}", test);
assert!(test.len() == 1);
assert_approx_eq!(test[0].0, 10.0);
assert_approx_eq!(test[0].1, 20.0);
}
} }

53
src/lib.rs
View File

@@ -47,16 +47,20 @@ pub struct Config {
pub start: f64, pub start: f64,
pub end: f64, pub end: f64,
pub uncorr: bool, pub uncorr: bool,
pub convdt: f64,
pub ignore_empty: bool
} }
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={:?}, write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?},
verbose={}, tolerance={}, iterations={}, temperature={}, verbose={}, tolerance={}, iterations={}, temperature={},
cyclic={:?}, uncorr={:?}, 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.metadata_file, self.hist_min, self.hist_max, self.num_bins,
self.verbose, self.tolerance, self.max_iterations, self.temperature, self.verbose, self.tolerance, self.max_iterations, self.temperature,
self.cyclic, self.uncorr, 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)
} }
} }
@@ -183,26 +187,39 @@ pub fn run(cfg: &Config) -> Result<()>{
println!("Supplied WHAM options: {}", &cfg); println!("Supplied WHAM options: {}", &cfg);
println!("Reading input files."); println!("Reading input files.");
let dataset = io::read_data(&cfg).chain_err(|| "Failed to create histogram.")?; let datasets = io::read_data(&cfg).chain_err(|| "Failed to read data.")?;
println!("{}", &dataset);
let (P, F, F_prev) = perform_wham(&cfg, &dataset)?; for (idx, dataset) in datasets.iter().enumerate() {
println!("WHAM converged."); 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 { let (P_std, free_energy_std) = if cfg.bootstrap > 0 {
println!("Bootstrapping.."); println!("Bootstrapping..");
error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap) error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap)
} else { } else {
(vec![0.0; P.len()], vec![0.0; P.len()]) (vec![0.0; P.len()], vec![0.0; P.len()])
}; };
// calculate free energy and dump state // calculate free energy and dump state
println!("Finished. Dumping final PMF"); println!("Finished. Dumping PMF");
let free_energy = calc_free_energy(&dataset, &P); 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) dump_state(&dataset, &F, &F_prev, &P, &P_std, &free_energy, &free_energy_std);
.chain_err(|| "Could not write results to output file")?; 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(()) Ok(())
} }

6
src/main.rs
View File

@@ -68,10 +68,14 @@ fn cli() -> Result<Config> {
} }
let dimens = num_bins.len(); 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, Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens,
verbose, tolerance, max_iterations, temperature, cyclic, output, verbose, tolerance, max_iterations, temperature, cyclic, output,
bootstrap, bootstrap_seed, start, end, uncorr}) bootstrap, bootstrap_seed, start, end, uncorr, convdt, ignore_empty})
} }
fn main() { fn main() {

108
tests/examples.rs
View File

@@ -6,64 +6,154 @@ mod integration {
use std::process::Command; use std::process::Command;
use std::fs; use std::fs;
use super::command::get_command; use super::command::get_command;
use std::fs::OpenOptions;
use std::io::prelude::*;
#[test] #[test]
fn wham_1d_cyclic() { fn wham_1d_cyclic() {
let output_file = "/tmp/wham_test_1d_cyclic.out";
get_command() get_command()
.args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"]) .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(&["-f", "example/1d_cyclic/metadata.dat"])
.args(&["-o", "/tmp/wham_test_1d_cyclic.out"]) .args(&["-o", output_file])
.output() .output()
.expect("failed to execute process"); .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") let output = Command::new("diff")
.arg("/tmp/wham_test_1d_cyclic.out") .arg(output_file)
.arg("example/1d_cyclic/wham.out") .arg("example/1d_cyclic/wham.out")
.output() .output()
.expect("failed to run diff"); .expect("failed to run diff");
let output_len = String::from_utf8_lossy(&output.stdout).len(); let output_len = String::from_utf8_lossy(&output.stdout).len();
assert_eq!(output_len, 0); assert_eq!(output_len, 0);
std::fs::remove_file(output_file).unwrap();
} }
#[test] #[test]
fn wham_1d_cyclic_uncorrelated() { fn wham_1d_cyclic_uncorrelated() {
let output_file = "/tmp/wham_test_1d_cyclic.out";
get_command() get_command()
.args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic", "--uncorr"]) .args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic", "--uncorr"])
.args(&["--seed", "1234"]) .args(&["--seed", "1234"])
.args(&["-f", "example/1d_cyclic/metadata.dat"]) .args(&["-f", "example/1d_cyclic/metadata.dat"])
.args(&["-o", "/tmp/wham_test_1d_cyclic.out"]) .args(&["-o", output_file])
.output() .output()
.expect("failed to execute process"); .expect("failed to execute process");
assert!(fs::metadata("/tmp/wham_test_1d_cyclic.out").is_ok()); assert!(fs::metadata("/tmp/wham_test_1d_cyclic.out").is_ok());
let output = Command::new("diff") let output = Command::new("diff")
.arg("/tmp/wham_test_1d_cyclic.out") .arg(output_file)
.arg("example/1d_cyclic/wham_uncorrelated.out") .arg("example/1d_cyclic/wham_uncorrelated.out")
.output() .output()
.expect("failed to run diff"); .expect("failed to run diff");
let output_len = String::from_utf8_lossy(&output.stdout).len(); let output_len = String::from_utf8_lossy(&output.stdout).len();
assert_eq!(output_len, 0); 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] #[test]
#[ignore] // expensive #[ignore] // expensive
fn wham_2d_cyclic() { fn wham_2d_cyclic() {
let output_file = "/tmp/wham_test_2d_cyclic.out";
get_command() get_command()
.args(&["--bins", "100,100", "--max", "pi,pi", "--min", "-pi,-pi", "-T", "300", "--cyclic"]) .args(&["--bins", "100,100", "--max", "pi,pi", "--min", "-pi,-pi", "-T", "300", "--cyclic"])
.args(&["-f", "example/2d_cyclic/metadata.dat"]) .args(&["-f", "example/2d_cyclic/metadata.dat"])
.args(&["-o", "/tmp/wham_test_2d_cyclic.out"]) .args(&["-o", output_file])
.output() .output()
.expect("failed to execute process"); .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") let output = Command::new("diff")
.arg("/tmp/wham_test_2d_cyclic.out") .arg(output_file)
.arg("example/2d_cyclic/wham.out") .arg("example/2d_cyclic/wham.out")
.output() .output()
.expect("failed to run diff"); .expect("failed to run diff");
let output_len = String::from_utf8_lossy(&output.stdout).len(); let output_len = String::from_utf8_lossy(&output.stdout).len();
assert_eq!(output_len, 0); assert_eq!(output_len, 0);
std::fs::remove_file(output_file).unwrap();
} }
} }