Update fep with alchemlyb and multiple edr files input #132
Conversation
mdpow/fep.py
Outdated
| @@ -693,7 +694,9 @@ def dgdl_edr(self, *args): | |||
| if not os.path.exists(fn): | |||
| logger.error("Missing dgdl.edr file %(fn)r.", vars()) | |||
| raise IOError(errno.ENOENT, "Missing dgdl.edr file", fn) | |||
| return fn | |||
| edrs = glob.glob(*args + (self.deffnm + '*.edr',)) | |||
There was a problem hiding this comment.
This looks as if it might wrongly raise an exception if there isn't a deffnm.edr file present.
How about something like (removing lines start ing at fn = ...)
pattern = os.path.join(args) + self.deffnm + '*.edr'
edrs = glob.glob(pattern)
if not edrs:
logger.error("Missing dgdl.edr file %(pattern)r.", vars())
raise IOError(errno.ENOENT, "Missing dgdl.edr file", pattern)
return [os.path.abspath(i) for i in edrs](Btw, I would be surprised if edrs = glob.glob(*args + (self.deffnm + '*.edr',)) actually worked.)
mdpow/fep.py
Outdated
| logger.info(" {0} --> {1}".format(edr, xvgfile)) | ||
| gromacs.g_energy(s=tpr, f=edr, odh=xvgfile) | ||
| dirct = os.path.abspath(os.path.dirname(edr[0])) | ||
| total_edr = os.path.join(dirct, 'total.edr') |
There was a problem hiding this comment.
Where does total.edr come from? Why is there a hard-coded name??
At least make the name configurable as a kwarg
def dgdl_tpr(self, *args, **kwargs):
total_edr_name = kwargs.get("total_edr_name", "total.edr")
mdpow/fep.py
Outdated
| dirct = os.path.abspath(os.path.dirname(edr[0])) | ||
| total_edr = os.path.join(dirct, 'total.edr') | ||
| logger.info(" {0} --> {1}".format('edrs', total_edr)) | ||
| total_edr = gromacs.eneconv(f=edr, o=total_edr) |
There was a problem hiding this comment.
What is total_edr now? I don't think that GW will return a specific filename. Should probably be
| total_edr = gromacs.eneconv(f=edr, o=total_edr) | |
| gromacs.eneconv(f=edr, o=total_edr) |
fix unit problem of VDW and Coulomb components
There was a problem hiding this comment.
Looking pretty good, especially with the tests. Some smaller things – please ping me when done & I merge and then look at mdpow-fep.
My main question is how you set self.method and self.SI and how you decide to use alchemlyb. That's important for knowing how to pass these values from the cfg file/scripts to the python API. If you add a short explanation in the comments or docs (even better) then I can attempt mdpow-fep.
|
Can you please look at #133 next? I don’t have time today to do it. |
-Use alchemlyb to calculate free energies with mbar and TI (if replace current TI calculation with alchemlyb, we should also deprecate plot function in
fep.py)-Allow multiple edr files for each fep window