mdpow-* script needed for analyzing BAR runs

[orbeckst]

At the moment there is no script available to calculate free energies from BAR runs. This is essential. If we need to call g_bar then that's ok but it must be wrapped in a mdpow-* script and fit into the standard workflow.

What's a workaround for right now?

[orbeckst]

Workaround: run g_bar and collect the total line:

total  0.000 -  1.000,   DG -1.09 +/-  0.21

For instance,

cd FEP/cyclohexane/VDW
g_bar -f */*.xvg  -o -oi -oh | awk '/^total/ {printf("%f %f\n", $6, $8)}'

In GromacsWrapper:

import gromacs
import glob
# capture command output
rc, stdout, stderr = gromacs.g_bar(f=glob.glob("*/*.xvg"), o="bar.xvg", oi="barint.xvg", oh="histogram.xvg", stdout=False, stderr=False)
# find and parse line 'total  0.000 -  1.000,   DG -1.09 +/-  0.21'
total = [line.split() for line in stdout.split('\n') if line.startswith('total')][-1]
DG, errDG = [float(total[field]) for field in (5, 7)]
print(DG, errDG)
# (-1.09, 0.21)

[orbeckst]

With #11 we have a general script and we can focus on it but it still cannot handle BAR.

[orbeckst]

@ianmkenney , the following is a summary of what we just discussed:

• mdpow-solvation should have an option --method that takes either "BAR" or "TI" and then uses the appropriate method to calculate the energies.
• I had the idea that the default should be whatever the runinput file says but the structure of the mdpow-solvation script is such that it does not look for runinput files but instead can be run on many different MDPOW directories. Thus, don't bother looking for the runiput value for the method.
• If one type of analysis (say BAR) is requested but no data are available, it should exit with an appropriate error.
• For a starting point, use g_bar for the calculation and error analysis. We can do other things (such as pymbar) later.
• In the code, start at fep.Gsolv.analysis(). Make the current one calc_DG_TI() and make a new calc_DG_BAR(). analyze(method=["BAR" | "TI"]) should then choose the appropriate calc_DG_*() method. In principle, it would be good to calculate both if possible but for right now, just calculate one or the other.
• Note that self.results.xvg is a dict containing keys Coulomb and VDW. Each contains a tuple with the list of lambdas and the corresponding XVG instance for each output file, XVG.real_filename should be the path to the file and can be used. However, if the file is compressed (ends in bz2 or gz) then you probably need to uncompress it (or does g_bar work with compressed files?)---see fep.Gsolv.compress_dgdl_xvg() for hints on how to do compression in Python (instead of system() calls...).

[orbeckst]

@iorga : note that this is issue is not a show-stopper because we save both TI and BAR data in the Gromacs run and we can process the TI data as before... at least I think that's we have in the xvg files: column 0 is time, column 1 is dH/dlambda (which is the one we parse for TI) and subsequent columns are the values of H for the neighboring lambdas.

From FEP/water/Coulomb/0500/md.xvg.bz2:

@    title "dH/d\xl\f{} and \xD\f{}H"
@    xaxis  label "Time (ps)"
@    yaxis  label "dH/d\xl\f{} and \xD\f{}H (kJ/mol [\xl\f{}]\S-1\N)"
@TYPE xy
@ subtitle "T = 300 (K) \xl\f{} = 0.5000"
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "dH/d\xl\f{} \xl\f{} 0.5000"
@ s1 legend "\xD\f{}H \xl\f{} 0.2500"
@ s2 legend "\xD\f{}H \xl\f{} 0.7500"
0.0000 155.55392 -38.888482 38.888482
0.1000 111.49819 -27.874549 27.874549

[orbeckst]

With PR #142 merged, alchemlyb MBAR is fully supported and the default.