[WIP] Issue 842 wham wrapper

[fiona-naughton]

WIP for #842

Uses the Grossfield WHAM implementation. Built on top of the auxiliary and bundle stuff in #868 and #900; only 'wham.py' is new here.

Given a datreant bundle (containing simulations from e.g. umbrella sampling), will create the necessary input files, run the Grossfield wham, and return the resulting profile (and remove the files).

The required data for each window should be stored in the bundle in categories (spring constant, value restrained to, ...) or as auxiliaries (the pull force/reaction coord data) under the same name for all simulations (the default name or a specified different name).

profile = wham(bundle)  # will look for default names
profile = wham(bundle, temperature='T', spring='k') # will look for provided names

The various command-line wham options (number of bins, etc) have been given default values (or in the case of the min/max reaction coordinate values, will be calculated) but can also be specified e.g. wham(bundle, num_bins=500, hist_min=1, hist_max=5). A start and end time can also be specified: wham(bundle, start_time=1000); only the datapoints within this time will be passed on to wham. Bootstrap error analysis is run by default ('profile' is a (num_bins, 3) array with reaction coord, PMF, and error) but can be switched off (w/ run_bootstrap=False; 'profile' is a (num_bins, 2) array); and the created files can be kept with keep_files.

Currently works for my test system; I don't have a test system for the different-temperature case but the input files appear to be generated correctly.

There's still a bit to do - particularly, getting the auxiliary data through to the input files feels like it could be tidier (I might play around with the auxiliary stuff to try make something nicer). A lot of the naming is to be in line with the Grossfield wham, but can probably be made clearer. Also need to deal with the fact the energy units may differ (Grossfield wham can be complied for kcal or kJ), but not sure what would be the best way to do that?

PR Checklist

• Tests?
• Docs?
• CHANGELOG updated?
• Issue raised/referenced?

[coveralls]

Coverage increased (+0.2%) to 84.354% when pulling 41f322d on fiona-naughton:issue-842-wham-wrapper into a1aa844 on MDAnalysis:develop.

[jbarnoud]

Glad to see this PR. It makes me realise that I really should merge #868. In the mean time, I can already tell you that your wham function is way too long and should be split.

[khuston]

FYI I tried to implement wham as described by Zhu and Hummer (2012), in which they avoid the slow iterative solution method in favor of using the BFGS algorithm to maximize the likelihood (minimize the negative log likelihood). A key advantage is supposed to be performance, cutting wham time by ~100x. Unfortunately the solver is going to crazy "solutions" that blow up the likelihood.

I haven't checked the behavior when I plug the correct solution as a guess into the solver, that would probably be the next step. The Jacobian function dAdg I copied from the paper gives the same Jacobian as scipy's minimizer gets from finite differences, so that's something. I can't spend any more time on it though, so just in case you'd like to pick it up and incorporate it. Obviously it's not much use in its current broken state:
https://github.com/khuston/pywham/blob/master/pywham.py

[coveralls]

Coverage increased (+0.2%) to 84.31% when pulling 5f94d20 on fiona-naughton:issue-842-wham-wrapper into 89f094e on MDAnalysis:develop.

[jbarnoud]

Couldn't this function be more generic and take a list of names as argument? Maybe something in the line of:

def check_names(bundle, mandatories, optionals)

[orbeckst]

@fiona-naughton how likely are you to pick up this work (and PR #900)?

One possibility would be to make it a small standalone module similar to MDAnalysis/RotamerConvolveMD. You would then have much more freedom to import MDS and the merge/review process would not be as stringent as for the MDAnalysis library itself. You can also decide how much tests you want.

This might be a way to get your existing code in a working state without tying up too much time. It would be a shame if all of your remaining work just languished as "PR: WIP".

[fiona-naughton]

@orbeckst Yes, I do keep meaning to get around to working more on this (if at least because I'm planning to use it myself!) - unfortunately in trying to catch up on other work it kept getting pushed back... When I get a chance to refresh my mind on what I was going I'll get back on this, hopefully in the near future! Thanks for the suggestion - that does sound like it would help

[orbeckst]

@fiona-naughton it looks that you have been busy with other things, which is totally fine. As I said above #923 (comment) and in PR #900, it would be great to have the wham functionality as a standalone module that also uses MDSynthesis. I am therefore closing this issue.

[fiona-naughton]

@orbeckst yes - I've always been meaning to get back to this but there always seemed to be other little things I wanted to finish first, sorry! As you say, it would still be nice as a standalone module; I did end up cobbling together some stuff based on this for my own analysis, and it'd be nice to talk more about it in my thesis, so hopefully that might finally motivate me to get something more tidy + usable.