From 0ebfb59def21c0f71afa9f462a1404919039cb3b Mon Sep 17 00:00:00 2001 From: VOD555 Date: Sat, 1 May 2021 12:54:46 -0700 Subject: [PATCH 01/20] add default setting for SI option --- scripts/mdpow-pcw | 2 +- scripts/mdpow-pow | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/mdpow-pcw b/scripts/mdpow-pcw index f760af91..ef89afe3 100755 --- a/scripts/mdpow-pcw +++ b/scripts/mdpow-pcw @@ -224,7 +224,7 @@ if __name__ == "__main__": "WRONG RESULTS CAN OCCUR! USE AT YOUR OWN RISK [%default]") parser.set_defaults(plotfile="none", outfile="pcw.txt", energyfile="energies.txt", - force=False, stride=1, permissive=False) + force=False, stride=1, permissive=False, SI=True) opts,args = parser.parse_args() if len(args) == 0: diff --git a/scripts/mdpow-pow b/scripts/mdpow-pow index a6d8c25f..f00b0dcb 100755 --- a/scripts/mdpow-pow +++ b/scripts/mdpow-pow @@ -224,7 +224,7 @@ if __name__ == "__main__": "WRONG RESULTS CAN OCCUR! USE AT YOUR OWN RISK [%default]") parser.set_defaults(plotfile="none", outfile="pow.txt", energyfile="energies.txt", - force=False, stride=1, permissive=False) + force=False, stride=1, permissive=False, SI=True) opts,args = parser.parse_args() if len(args) == 0: From 092b9838189cd0f7afbc77a0b6bc4da459799594 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Sat, 1 May 2021 15:28:06 -0700 Subject: [PATCH 02/20] update docs --- scripts/mdpow-pcw | 8 ++++++-- scripts/mdpow-pow | 8 ++++++-- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/scripts/mdpow-pcw b/scripts/mdpow-pcw index ef89afe3..c03560b9 100755 --- a/scripts/mdpow-pcw +++ b/scripts/mdpow-pcw @@ -190,7 +190,9 @@ if __name__ == "__main__": parser.add_option('-p', '--plotfile', dest="plotfile", help="plot dV/dlambda to FILE; use png or pdf suffix to " "determine the file type. 'auto' generates a pdf file " - "DIRECTORY/dVdl_MOLID_pcw.pdf and 'None' disables it [%default]", + "DIRECTORY/dVdl_MOLID_pcw.pdf and 'None' disables it [%default]" + "The plot function is only available for mdpow estimator," + "and is disabled when using alchemlyb estimators.", metavar="FILE") parser.add_option('-o', '--outfile', dest="outfile", help="append one-line results summary to FILE [%default]", @@ -207,7 +209,9 @@ if __name__ == "__main__": help="force rereading all data [%default]") parser.add_option('--noSI', dest='SI', action="store_false", - help="Apply statistical inefficiency analysis") + help="Disable statistical inefficiency analysis." + "Statitical inefficiency analysis is performed by default when using" + "alchemlyb estimators and is disabled when using mdpow estimator.") parser.add_option('-s', '--stride', dest="stride", type='int', help="Use every N-th datapoint from the original dV/dlambda data. " "[%default]", diff --git a/scripts/mdpow-pow b/scripts/mdpow-pow index f00b0dcb..86649ca6 100755 --- a/scripts/mdpow-pow +++ b/scripts/mdpow-pow @@ -190,7 +190,9 @@ if __name__ == "__main__": parser.add_option('-p', '--plotfile', dest="plotfile", help="plot dV/dlambda to FILE; use png or pdf suffix to " "determine the file type. 'auto' generates a pdf file " - "DIRECTORY/dVdl_MOLID_pow.pdf and 'None' disables it [%default]", + "DIRECTORY/dVdl_MOLID_pow.pdf and 'None' disables it [%default]" + "The plot function is only available for mdpow estimator," + "and is disabled when using alchemlyb estimators.", metavar="FILE") parser.add_option('-o', '--outfile', dest="outfile", help="append one-line results summary to FILE [%default]", @@ -207,7 +209,9 @@ if __name__ == "__main__": help="force rereading all data [%default]") parser.add_option('--noSI', dest='SI', action="store_false", - help="Apply statistical inefficiency analysis") + help="Disable statistical inefficiency analysis" + "Statitical inefficiency analysis is performed by default when using" + "alchemlyb estimators and is disabled when using mdpow estimator.") parser.add_option('-s', '--stride', dest="stride", type='int', help="Use every N-th datapoint from the original dV/dlambda data. " "[%default]", From f74ff8236a609ca70befb70951efbbeec3c8f68c Mon Sep 17 00:00:00 2001 From: VOD555 Date: Wed, 2 Jun 2021 07:16:47 -0700 Subject: [PATCH 03/20] add new options for solvationenergy --- mdpow/fep.py | 8 ++-- scripts/mdpow-solvationenergy | 86 +++++++++++++++++++++++++++++------ 2 files changed, 77 insertions(+), 17 deletions(-) diff --git a/mdpow/fep.py b/mdpow/fep.py index c833f055..793fa157 100644 --- a/mdpow/fep.py +++ b/mdpow/fep.py @@ -1368,14 +1368,14 @@ def p_transfer(G1, G2, **kwargs): G.method = kwargs.pop('method', 'MBAR') if kwargs['force']: if estimator == 'mdpow': - G.start = kwargs.pop('start', 0) - G.stop = kwargs.pop('stop', None) - G.SI = kwargs.pop('stop', False) - G.analyze(**kwargs) if G.method != "TI": errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method logger.error(errmsg) raise ValueError(errmsg) + G.start = kwargs.pop('start', 0) + G.stop = kwargs.pop('stop', None) + G.SI = kwargs.pop('stop', False) + G.analyze(**kwargs) elif estimator == 'alchemlyb': G.analyze_alchemlyb(**kwargs) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index 13d9dd36..e922e78e 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -105,11 +105,39 @@ def run_gsolv(solvent, directory, **kwargs): gsolv = load_gsolv(solvent, directory, permissive=kwargs.pop('permissive',False),) + if not hasattr(gsolv, 'start'): + gsolv.start = kwargs.pop('start', 0) + if not hasattr(gsolv, 'stop'): + gsolv.stop = kwargs.pop('stop', None) + if not hasattr(gsolv, 'SI'): + kwargs.setdefault('SI', True) + else: + kwargs.setdefault('SI', gsolv.SI) + + # for this version. use the method given instead of the one in the input cfg file + gsolv.method = kwargs.pop('method', 'MBAR') + estimator = kwargs.pop('estimator', 'alchemlyb') + # do the analysis (only use keywords that are needed for FEP.analysis()) akw = {'stride': kwargs.pop('stride', None), 'force': kwargs.pop('force', False)} + if not estimator in ('mdpow', 'alchemlyb'): + errmsg = "estimator = %r is not supported, must be 'mdpow' or 'alchemlyb'" % estimator + logger.error(errmsg) + raise ValueError(errmsg) if akw['force']: logger.info("[%(directory)s] Forcing (re)analysis of hydration free energy data", vars()) - gsolv.analyze(**akw) + if estimator == 'mdpow': + if G.method != "TI": + errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method + logger.error(errmsg) + raise ValueError(errmsg) + gsolv.start = kwargs.pop('start', 0) + gsolv.stop = kwargs.pop('stop', None) + gsolv.SI = kwargs.pop('stop', False) + gsolv.analyze(**akw) + + elif estimator == 'alchemlyb': + gsolv.analyze_alchemlyb(**akw) # make sure that we have data try: gsolv.results.DeltaA.Gibbs @@ -117,8 +145,19 @@ def run_gsolv(solvent, directory, **kwargs): gsolv.logger_DeltaA0() except (KeyError, AttributeError): # KeyError because results is a AttributeDict logger.info("[%(directory)s] Analyzing %(solvent)s solvation free energy data for the first time", vars()) - gsolv.analyze(**akw) - + if estimator == 'mdpow': + if G.method != "TI": + errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method + logger.error(errmsg) + raise ValueError(errmsg) + gsolv.start = kwargs.pop('start', 0) + gsolv.stop = kwargs.pop('stop', None) + gsolv.SI = kwargs.pop('stop', False) + gsolv.analyze(**akw) + + elif estimator == 'alchemlyb': + gsolv.analyze_alchemlyb(**akw) + if datstat(gsolv) == 'BAD': logger.warning("[%s] Possible file corruption: %s:%s", directory, solvent, datstat(gsolv)) @@ -129,6 +168,8 @@ def run_gsolv(solvent, directory, **kwargs): logger.info("[%s] Wrote solvation energy terms to %r", directory, kwargs['energyfile']) plotfile = kwargs.get('plotfile', None) + if estimator == 'alchemlyb': + plotfile = False if plotfile: if plotfile == 'auto': plotfile = auto_plotfile(directory, gsolv) @@ -185,31 +226,48 @@ if __name__ == "__main__": parser.add_option('-p', '--plotfile', dest="plotfile", help="plot dV/dlambda to FILE; use png or pdf suffix to " "determine the file type. 'auto' generates a pdf file " - "DIRECTORY/dVdl__.pdf and 'None' disables it [%default]", + "DIRECTORY/dVdl_MOLID_pow.pdf and 'None' disables it [%default]" + "The plot function is only available for mdpow estimator," + "and is disabled when using alchemlyb estimators.", + metavar="FILE") + parser.add_option('--solvent', '-S', dest='solvent', + metavar="NAME", + help="solvent NAME for compound, 'water', 'octanol', or 'cyclohexane' [%default]") + parser.add_option('-o', '--outfile', dest="outfile", + help="append one-line results summary to FILE [%default]", metavar="FILE") parser.add_option('-e', '--energies', dest="energyfile", help="append solvation free energies to FILE [%default]", metavar="FILE") + parser.add_option('--estimator', dest="estimator", default='alchemlyb', + help="choose the estimator to be used, alchemlyb or mdpow estimators") + parser.add_option('--method', dest="method", default='MBAR', + help="choose the method to calculate free energy",) parser.add_option('--force', dest='force', action="store_true", help="force rereading all data [%default]") + parser.add_option('--noSI', dest='SI', + action="store_false", + help="Disable statistical inefficiency analysis" + "Statitical inefficiency analysis is performed by default when using" + "alchemlyb estimators and is disabled when using mdpow estimator.") parser.add_option('-s', '--stride', dest="stride", type='int', help="Use every N-th datapoint from the original dV/dlambda data. " - "Using a number larger than 1 can significantly reduce memory usage " - "with little impact on the precision of the final output. [%default]", + "[%default]", metavar="N") - parser.add_option('--solvent', '-S', dest='solvent', - metavar="NAME", - help="solvent NAME for compound, 'water', 'octanol', or 'cyclohexane' [%default]") + parser.add_option('--start', dest="start", type='int', + help="Start point for the data used from the original dV/dlambda data. ") + parser.add_option('--stop', dest="stop", type='int', + help="Stop point for the data used from the original dV/dlambda data.",) parser.add_option('--ignore-corrupted', dest="permissive", action="store_true", help="skip lines in the md.xvg files that cannot be parsed. " "WARNING: Other lines in the file might have been corrupted in " "such a way that they appear correct but are in fact wrong. " "WRONG RESULTS CAN OCCUR! USE AT YOUR OWN RISK [%default]") - parser.set_defaults(plotfile="auto", - energyfile="gsolv.txt", - force=False, stride=10, permissive=False, solvent='water') + parser.set_defaults(plotfile="none", + outfile="pow.txt", energyfile="energies.txt", + force=False, stride=1, permissive=False, SI=True) opts,args = parser.parse_args() if len(args) == 0: @@ -230,7 +288,9 @@ if __name__ == "__main__": try: run_gsolv(opts.solvent, directory, plotfile=opts.plotfile, energyfile=realpath_or_none(opts.energyfile), - force=opts.force, stride=opts.stride, permissive=opts.permissive) + force=opts.force, stride=opts.stride, start=opts.start, + stop=opts.stop, estimator=opts.estimator, + method=opts.method, permissive=opts.permissive, SI=opts.SI) except (OSError, IOError), err: logger.error("Running analysis in directory %r failed: %s", directory, str(err)) except Exception, err: From 33d1511aa6e078d59162989ddb125ac699f3e501 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Wed, 2 Jun 2021 20:14:07 -0700 Subject: [PATCH 04/20] use argparse instead of optionparser --- scripts/mdpow-solvationenergy | 145 +++++++++++++++------------------- 1 file changed, 65 insertions(+), 80 deletions(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index e922e78e..122540f9 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -105,6 +105,19 @@ def run_gsolv(solvent, directory, **kwargs): gsolv = load_gsolv(solvent, directory, permissive=kwargs.pop('permissive',False),) + # for this version. use the method given instead of the one in the input cfg file + gsolv.method = kwargs.pop('method', 'MBAR') + estimator = kwargs.pop('estimator', 'alchemlyb') + if not estimator in ('mdpow', 'alchemlyb'): + errmsg = "estimator = %r is not supported, must be 'mdpow' or 'alchemlyb'" % estimator + logger.error(errmsg) + raise ValueError(errmsg) + if estimator == 'mdpow': + if gsolv.method != "TI": + errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method + logger.error(errmsg) + raise ValueError(errmsg) + if not hasattr(gsolv, 'start'): gsolv.start = kwargs.pop('start', 0) if not hasattr(gsolv, 'stop'): @@ -114,50 +127,20 @@ def run_gsolv(solvent, directory, **kwargs): else: kwargs.setdefault('SI', gsolv.SI) - # for this version. use the method given instead of the one in the input cfg file - gsolv.method = kwargs.pop('method', 'MBAR') - estimator = kwargs.pop('estimator', 'alchemlyb') - # do the analysis (only use keywords that are needed for FEP.analysis()) akw = {'stride': kwargs.pop('stride', None), 'force': kwargs.pop('force', False)} - if not estimator in ('mdpow', 'alchemlyb'): - errmsg = "estimator = %r is not supported, must be 'mdpow' or 'alchemlyb'" % estimator - logger.error(errmsg) - raise ValueError(errmsg) - if akw['force']: - logger.info("[%(directory)s] Forcing (re)analysis of hydration free energy data", vars()) - if estimator == 'mdpow': - if G.method != "TI": - errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method - logger.error(errmsg) - raise ValueError(errmsg) - gsolv.start = kwargs.pop('start', 0) - gsolv.stop = kwargs.pop('stop', None) - gsolv.SI = kwargs.pop('stop', False) - gsolv.analyze(**akw) - elif estimator == 'alchemlyb': - gsolv.analyze_alchemlyb(**akw) # make sure that we have data - try: - gsolv.results.DeltaA.Gibbs + if akw['force'] or (not hasattr(gsolv.results.DeltaA, 'Gibbs')): + logger.info("[%(directory)s] Forcing (re)analysis of hydration free energy data", vars()) logger.info("[%(directory)s] Reading existing data from pickle file.", vars()) - gsolv.logger_DeltaA0() - except (KeyError, AttributeError): # KeyError because results is a AttributeDict - logger.info("[%(directory)s] Analyzing %(solvent)s solvation free energy data for the first time", vars()) if estimator == 'mdpow': - if G.method != "TI": - errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method - logger.error(errmsg) - raise ValueError(errmsg) - gsolv.start = kwargs.pop('start', 0) - gsolv.stop = kwargs.pop('stop', None) - gsolv.SI = kwargs.pop('stop', False) gsolv.analyze(**akw) - elif estimator == 'alchemlyb': gsolv.analyze_alchemlyb(**akw) - + + gsolv.logger_DeltaA0() + if datstat(gsolv) == 'BAD': logger.warning("[%s] Possible file corruption: %s:%s", directory, solvent, datstat(gsolv)) @@ -220,51 +203,53 @@ def realpath_or_none(fn): if __name__ == "__main__": import sys import os.path - from optparse import OptionParser - - parser = OptionParser(usage=__doc__) - parser.add_option('-p', '--plotfile', dest="plotfile", - help="plot dV/dlambda to FILE; use png or pdf suffix to " - "determine the file type. 'auto' generates a pdf file " - "DIRECTORY/dVdl_MOLID_pow.pdf and 'None' disables it [%default]" - "The plot function is only available for mdpow estimator," - "and is disabled when using alchemlyb estimators.", - metavar="FILE") - parser.add_option('--solvent', '-S', dest='solvent', - metavar="NAME", - help="solvent NAME for compound, 'water', 'octanol', or 'cyclohexane' [%default]") - parser.add_option('-o', '--outfile', dest="outfile", - help="append one-line results summary to FILE [%default]", - metavar="FILE") - parser.add_option('-e', '--energies', dest="energyfile", - help="append solvation free energies to FILE [%default]", - metavar="FILE") - parser.add_option('--estimator', dest="estimator", default='alchemlyb', - help="choose the estimator to be used, alchemlyb or mdpow estimators") - parser.add_option('--method', dest="method", default='MBAR', - help="choose the method to calculate free energy",) - parser.add_option('--force', dest='force', - action="store_true", - help="force rereading all data [%default]") - parser.add_option('--noSI', dest='SI', - action="store_false", - help="Disable statistical inefficiency analysis" - "Statitical inefficiency analysis is performed by default when using" - "alchemlyb estimators and is disabled when using mdpow estimator.") - parser.add_option('-s', '--stride', dest="stride", type='int', - help="Use every N-th datapoint from the original dV/dlambda data. " - "[%default]", - metavar="N") - parser.add_option('--start', dest="start", type='int', - help="Start point for the data used from the original dV/dlambda data. ") - parser.add_option('--stop', dest="stop", type='int', - help="Stop point for the data used from the original dV/dlambda data.",) - parser.add_option('--ignore-corrupted', dest="permissive", - action="store_true", - help="skip lines in the md.xvg files that cannot be parsed. " - "WARNING: Other lines in the file might have been corrupted in " - "such a way that they appear correct but are in fact wrong. " - "WRONG RESULTS CAN OCCUR! USE AT YOUR OWN RISK [%default]") + import argparse + + parser = argparse.ArgumentParser(usage=__doc__) + parser.add_argument('-p', '--plotfile', dest="plotfile", + help="plot dV/dlambda to FILE; use png or pdf suffix to " + "determine the file type. 'auto' generates a pdf file " + "DIRECTORY/dVdl__.pdf and 'None' disables it [%default]" + "The plot function is only available for mdpow estimator," + "and is disabled when using alchemlyb estimators.", + metavar="FILE") + parser.add_argument('--solvent', '-S', dest='solvent', + choices=['water', 'octanol', 'cyclohexane'], metavar="NAME", + help="solvent NAME for compound, 'water', 'octanol', or 'cyclohexane' [%default]") + parser.add_argument('-o', '--outfile', dest="outfile", + help="append one-line results summary to FILE [%default]", + metavar="FILE") + parser.add_argument('-e', '--energies', dest="energyfile", + help="append solvation free energies to FILE [%default]", + metavar="FILE") + parser.add_argument('--estimator', dest="estimator", default='alchemlyb', + choices=['mdpow', 'alchemlyb'], + help="choose the estimator to be used, alchemlyb or mdpow estimators") + parser.add_argument('--method', dest="method", default='MBAR', + choices=['TI', 'MBAR', 'BAR'], + help="choose the method to calculate free energy",) + parser.add_argument('--force', dest='force', + action="store_true", + help="force rereading all data [%default]") + parser.add_argument('--noSI', dest='SI', + action="store_false", + help="Disable statistical inefficiency analysis" + "Statitical inefficiency analysis is performed by default when using" + "alchemlyb estimators and is disabled when using mdpow estimator.") + parser.add_argument('-s', '--stride', dest="stride", type='int', + help="Use every N-th datapoint from the original dV/dlambda data. " + "[%default]", + metavar="N") + parser.add_argument('--start', dest="start", type='int', + help="Start point for the data used from the original dV/dlambda data. ") + parser.add_argument('--stop', dest="stop", type='int', + help="Stop point for the data used from the original dV/dlambda data.",) + parser.add_argument('--ignore-corrupted', dest="permissive", + action="store_true", + help="skip lines in the md.xvg files that cannot be parsed. " + "WARNING: Other lines in the file might have been corrupted in " + "such a way that they appear correct but are in fact wrong. " + "WRONG RESULTS CAN OCCUR! USE AT YOUR OWN RISK [%default]") parser.set_defaults(plotfile="none", outfile="pow.txt", energyfile="energies.txt", force=False, stride=1, permissive=False, SI=True) From 3851606da7a42dc399b8e1653d8302e15db611c1 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Wed, 2 Jun 2021 20:14:43 -0700 Subject: [PATCH 05/20] remove duplicate initial settings --- mdpow/fep.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/mdpow/fep.py b/mdpow/fep.py index 793fa157..e56fda70 100644 --- a/mdpow/fep.py +++ b/mdpow/fep.py @@ -1372,9 +1372,6 @@ def p_transfer(G1, G2, **kwargs): errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method logger.error(errmsg) raise ValueError(errmsg) - G.start = kwargs.pop('start', 0) - G.stop = kwargs.pop('stop', None) - G.SI = kwargs.pop('stop', False) G.analyze(**kwargs) elif estimator == 'alchemlyb': G.analyze_alchemlyb(**kwargs) From 905aa5e6b77bc8ee05e2c352823ef0fa276f2fd2 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 3 Jun 2021 05:58:34 -0700 Subject: [PATCH 06/20] use positional args for directory input --- scripts/mdpow-solvationenergy | 60 +++++++++++++++++------------------ 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index 122540f9..6883887e 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -1,5 +1,5 @@ #!/usr/bin/env python -"""%prog [options] DIRECTORY [DIRECTORY ...] +"""%(prog)s [options] DIRECTORY [DIRECTORY ...] Run the free energy analysis for a solvent in /FEP and return DeltaG. @@ -205,64 +205,64 @@ if __name__ == "__main__": import os.path import argparse - parser = argparse.ArgumentParser(usage=__doc__) - parser.add_argument('-p', '--plotfile', dest="plotfile", + parser = argparse.ArgumentParser(usage=__doc__, prog='mdpow-solvationenergy', + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument("directory", nargs='+', + help="directory or directories which contain all the files " + "resulting from running mdpow.fep.Ghyd.setup() ",) + parser.add_argument('--plotfile', dest="plotfile", default='none', help="plot dV/dlambda to FILE; use png or pdf suffix to " "determine the file type. 'auto' generates a pdf file " - "DIRECTORY/dVdl__.pdf and 'None' disables it [%default]" + "DIRECTORY/dVdl__.pdf and 'None' disables it " "The plot function is only available for mdpow estimator," "and is disabled when using alchemlyb estimators.", metavar="FILE") - parser.add_argument('--solvent', '-S', dest='solvent', + parser.add_argument('--solvent', '-S', dest='solvent', default='water', choices=['water', 'octanol', 'cyclohexane'], metavar="NAME", - help="solvent NAME for compound, 'water', 'octanol', or 'cyclohexane' [%default]") - parser.add_argument('-o', '--outfile', dest="outfile", - help="append one-line results summary to FILE [%default]", + help="solvent NAME for compound, 'water', 'octanol', or 'cyclohexane' ") + parser.add_argument('-o', '--outfile', dest="outfile", default="gsolv.txt", + help="append one-line results summary to FILE ", metavar="FILE") - parser.add_argument('-e', '--energies', dest="energyfile", - help="append solvation free energies to FILE [%default]", + parser.add_argument('-e', '--energies', dest="energyfile", default="energies.txt", + help="append solvation free energies to FILE ", metavar="FILE") parser.add_argument('--estimator', dest="estimator", default='alchemlyb', choices=['mdpow', 'alchemlyb'], help="choose the estimator to be used, alchemlyb or mdpow estimators") parser.add_argument('--method', dest="method", default='MBAR', choices=['TI', 'MBAR', 'BAR'], - help="choose the method to calculate free energy",) - parser.add_argument('--force', dest='force', + help="choose the method to calculate free energy ",) + parser.add_argument('--force', dest='force', default=False, action="store_true", - help="force rereading all data [%default]") + help="force rereading all data ") parser.add_argument('--noSI', dest='SI', - action="store_false", + action="store_false", default=True, help="Disable statistical inefficiency analysis" "Statitical inefficiency analysis is performed by default when using" - "alchemlyb estimators and is disabled when using mdpow estimator.") - parser.add_argument('-s', '--stride', dest="stride", type='int', - help="Use every N-th datapoint from the original dV/dlambda data. " - "[%default]", + "alchemlyb estimators and is disabled when using mdpow estimator. ") + parser.add_argument('-s', '--stride', dest="stride", type=int, default=1, + help="Use every N-th datapoint from the original dV/dlambda data. ", metavar="N") - parser.add_argument('--start', dest="start", type='int', + parser.add_argument('--start', dest="start", type=int, help="Start point for the data used from the original dV/dlambda data. ") - parser.add_argument('--stop', dest="stop", type='int', - help="Stop point for the data used from the original dV/dlambda data.",) + parser.add_argument('--stop', dest="stop", type=int, + help="Stop point for the data used from the original dV/dlambda data. ",) parser.add_argument('--ignore-corrupted', dest="permissive", - action="store_true", + action="store_true", default=False, help="skip lines in the md.xvg files that cannot be parsed. " "WARNING: Other lines in the file might have been corrupted in " "such a way that they appear correct but are in fact wrong. " - "WRONG RESULTS CAN OCCUR! USE AT YOUR OWN RISK [%default]") - parser.set_defaults(plotfile="none", - outfile="pow.txt", energyfile="energies.txt", - force=False, stride=1, permissive=False, SI=True) - opts,args = parser.parse_args() + "WRONG RESULTS CAN OCCUR! USE AT YOUR OWN RISK ") + opts = parser.parse_args() - if len(args) == 0: + if len(opts.directory) == 0: logger.fatal("A directory is required. See --help.") sys.exit(1) - elif len(args) > 1 and not opts.plotfile.lower() in ('none', 'auto'): + elif len(opts.directory) > 1 and not opts.plotfile.lower() in ('none', 'auto'): logger.fatal("Can only use --plotfile=None or --plotfile=auto with multiple directories.") sys.exit(1) - for directory in args: + for directory in opts.directory: if not os.path.exists(directory): logger.warn("Directory %r not found, skipping...", directory) continue From a8e3787d48400aefe0535c4886fc0ef36b34b16d Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 3 Jun 2021 11:45:26 -0700 Subject: [PATCH 07/20] make opts.plotfile unchanged during the loop for multiple directory inputs --- scripts/mdpow-solvationenergy | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index 6883887e..6aac2146 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -268,10 +268,10 @@ if __name__ == "__main__": continue logger.info("Analyzing directory %r... (can take a while)", directory) - opts.plotfile = parse_filename_option(opts.plotfile) + plotfile = parse_filename_option(opts.plotfile) try: - run_gsolv(opts.solvent, directory, plotfile=opts.plotfile, + run_gsolv(opts.solvent, directory, plotfile=plotfile, energyfile=realpath_or_none(opts.energyfile), force=opts.force, stride=opts.stride, start=opts.start, stop=opts.stop, estimator=opts.estimator, From 1dba7da402e893cc09a6ac69b68a575954b3310d Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 3 Jun 2021 21:01:37 -0700 Subject: [PATCH 08/20] add INFO loggers --- scripts/mdpow-solvationenergy | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index 6aac2146..cbe8ef13 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -104,17 +104,19 @@ def run_gsolv(solvent, directory, **kwargs): return _datastatus[g.contains_corrupted_xvgs()] gsolv = load_gsolv(solvent, directory, permissive=kwargs.pop('permissive',False),) + logger.info("The solvent is %s .", solvent) # for this version. use the method given instead of the one in the input cfg file - gsolv.method = kwargs.pop('method', 'MBAR') estimator = kwargs.pop('estimator', 'alchemlyb') if not estimator in ('mdpow', 'alchemlyb'): errmsg = "estimator = %r is not supported, must be 'mdpow' or 'alchemlyb'" % estimator logger.error(errmsg) raise ValueError(errmsg) + + gsolv.method = kwargs.pop('method', 'MBAR') if estimator == 'mdpow': if gsolv.method != "TI": - errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method + errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % gsolv.method logger.error(errmsg) raise ValueError(errmsg) @@ -133,11 +135,22 @@ def run_gsolv(solvent, directory, **kwargs): # make sure that we have data if akw['force'] or (not hasattr(gsolv.results.DeltaA, 'Gibbs')): logger.info("[%(directory)s] Forcing (re)analysis of hydration free energy data", vars()) - logger.info("[%(directory)s] Reading existing data from pickle file.", vars()) + # write out the settings when the analysis is performed + logger.info("Estimator is %s .", estimator) + logger.info("Free energy calculation method is %s .", gsolv.method) + logger.info("Analysis starts from frame %s .", gsolv.start) + logger.info("Analysis stops at frame %s .", gsolv.stop) + logger.info("Analysis stride is %s .", akw['stride']) + if kwargs['SI']: + logger.info("Statistial inefficiency analysis will be performed.") + else: + logger.info("Statistial inefficiency analysis won't be performed.") if estimator == 'mdpow': gsolv.analyze(**akw) elif estimator == 'alchemlyb': gsolv.analyze_alchemlyb(**akw) + else: + logger.info("[%(directory)s] Reading existing data from pickle file.", vars()) gsolv.logger_DeltaA0() @@ -243,9 +256,9 @@ if __name__ == "__main__": parser.add_argument('-s', '--stride', dest="stride", type=int, default=1, help="Use every N-th datapoint from the original dV/dlambda data. ", metavar="N") - parser.add_argument('--start', dest="start", type=int, + parser.add_argument('--start', dest="start", type=int, default=0, help="Start point for the data used from the original dV/dlambda data. ") - parser.add_argument('--stop', dest="stop", type=int, + parser.add_argument('--stop', dest="stop", type=int, default=None, help="Stop point for the data used from the original dV/dlambda data. ",) parser.add_argument('--ignore-corrupted', dest="permissive", action="store_true", default=False, From d3f78d88c3e123f6be4f912cce68431a3e45fb59 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Fri, 4 Jun 2021 02:23:25 -0700 Subject: [PATCH 09/20] fix description for noSI --- scripts/mdpow-solvationenergy | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index cbe8ef13..4844f7db 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -136,11 +136,11 @@ def run_gsolv(solvent, directory, **kwargs): if akw['force'] or (not hasattr(gsolv.results.DeltaA, 'Gibbs')): logger.info("[%(directory)s] Forcing (re)analysis of hydration free energy data", vars()) # write out the settings when the analysis is performed - logger.info("Estimator is %s .", estimator) - logger.info("Free energy calculation method is %s .", gsolv.method) - logger.info("Analysis starts from frame %s .", gsolv.start) - logger.info("Analysis stops at frame %s .", gsolv.stop) - logger.info("Analysis stride is %s .", akw['stride']) + logger.info("Estimator is %s.", estimator) + logger.info("Free energy calculation method is %s.", gsolv.method) + logger.info("Analysis starts from frame %s.", gsolv.start) + logger.info("Analysis stops at frame %s.", gsolv.stop) + logger.info("Analysis stride is %s.", akw['stride']) if kwargs['SI']: logger.info("Statistial inefficiency analysis will be performed.") else: @@ -250,16 +250,16 @@ if __name__ == "__main__": help="force rereading all data ") parser.add_argument('--noSI', dest='SI', action="store_false", default=True, - help="Disable statistical inefficiency analysis" + help="disable statistical inefficiency analysis. " "Statitical inefficiency analysis is performed by default when using" "alchemlyb estimators and is disabled when using mdpow estimator. ") parser.add_argument('-s', '--stride', dest="stride", type=int, default=1, - help="Use every N-th datapoint from the original dV/dlambda data. ", + help="use every N-th datapoint from the original dV/dlambda data. ", metavar="N") parser.add_argument('--start', dest="start", type=int, default=0, - help="Start point for the data used from the original dV/dlambda data. ") + help="start point for the data used from the original dV/dlambda data. ") parser.add_argument('--stop', dest="stop", type=int, default=None, - help="Stop point for the data used from the original dV/dlambda data. ",) + help="stop point for the data used from the original dV/dlambda data. ",) parser.add_argument('--ignore-corrupted', dest="permissive", action="store_true", default=False, help="skip lines in the md.xvg files that cannot be parsed. " From ca2ef03cd547fb491b5c7080632d9b5fa076e74d Mon Sep 17 00:00:00 2001 From: VOD555 Date: Fri, 4 Jun 2021 02:23:48 -0700 Subject: [PATCH 10/20] print out lambda window info when doing SI analysis --- mdpow/fep.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mdpow/fep.py b/mdpow/fep.py index e56fda70..8acb677f 100644 --- a/mdpow/fep.py +++ b/mdpow/fep.py @@ -1021,7 +1021,7 @@ def collect_alchemlyb(self, SI=True, start=0, stop=None, stride=None, autosave=T xvg_file = self.dgdl_xvg(self.wdir(component, l)) xvg_df = extract(xvg_file, T=self.Temperature).iloc[start:stop:stride] if SI: - logger.info("Perform statistical inefficiency analysis.") + logger.info("Performing statistical inefficiency analysis for window %s %s.", component, l) ts = _extract_dataframe(xvg_file).iloc[start:stop:stride] ts = pd.DataFrame({'time': ts.iloc[:,0], 'dhdl': ts.iloc[:,1]}) ts = ts.set_index('time') From eabc7eae70a94fd1c11c80b0f82b9660921b2c00 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Fri, 4 Jun 2021 02:53:17 -0700 Subject: [PATCH 11/20] disable alchemlyb loggers when using mdpow estimator --- scripts/mdpow-solvationenergy | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index 4844f7db..6b13e5cd 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -138,16 +138,17 @@ def run_gsolv(solvent, directory, **kwargs): # write out the settings when the analysis is performed logger.info("Estimator is %s.", estimator) logger.info("Free energy calculation method is %s.", gsolv.method) - logger.info("Analysis starts from frame %s.", gsolv.start) - logger.info("Analysis stops at frame %s.", gsolv.stop) logger.info("Analysis stride is %s.", akw['stride']) - if kwargs['SI']: - logger.info("Statistial inefficiency analysis will be performed.") - else: - logger.info("Statistial inefficiency analysis won't be performed.") if estimator == 'mdpow': gsolv.analyze(**akw) elif estimator == 'alchemlyb': + # the options only available for alchemlyb + logger.info("Analysis starts from frame %s.", gsolv.start) + logger.info("Analysis stops at frame %s.", gsolv.stop) + if kwargs['SI']: + logger.info("Statistial inefficiency analysis will be performed.") + else: + logger.info("Statistial inefficiency analysis won't be performed.") gsolv.analyze_alchemlyb(**akw) else: logger.info("[%(directory)s] Reading existing data from pickle file.", vars()) From 3a62bcd796873756b888b0df6ae70c7edd65ea8e Mon Sep 17 00:00:00 2001 From: VOD555 Date: Fri, 4 Jun 2021 02:53:40 -0700 Subject: [PATCH 12/20] add loggers for calculation setting in fep.py --- mdpow/fep.py | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/mdpow/fep.py b/mdpow/fep.py index 8acb677f..5c5d3131 100644 --- a/mdpow/fep.py +++ b/mdpow/fep.py @@ -1346,6 +1346,7 @@ def p_transfer(G1, G2, **kwargs): errmsg = "estimator = %r is not supported, must be 'mdpow' or 'alchemlyb'" % estimator logger.error(errmsg) raise ValueError(errmsg) + if G1.molecule != G2.molecule: raise ValueError("The two simulations were done for different molecules.") if G1.Temperature != G2.Temperature: @@ -1366,14 +1367,26 @@ def p_transfer(G1, G2, **kwargs): # for this version. use the method given instead of the one in the input cfg file G.method = kwargs.pop('method', 'MBAR') - if kwargs['force']: + if estimator == 'mdpow': + if G.method != "TI": + errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method + logger.error(errmsg) + raise ValueError(errmsg) + + if kwargs['force'] or (not hasattr(G.results.DeltaA, 'Gibbs')): + # write out the settings when the analysis is performed + logger.info("Estimator is %s.", estimator) + logger.info("Free energy calculation method is %s.", G.method) + logger.info("Analysis stride is %s.", akw['stride']) if estimator == 'mdpow': - if G.method != "TI": - errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method - logger.error(errmsg) - raise ValueError(errmsg) G.analyze(**kwargs) elif estimator == 'alchemlyb': + logger.info("Analysis starts from frame %s.", G.start) + logger.info("Analysis stops at frame %s.", G.stop) + if kwargs['SI']: + logger.info("Statistial inefficiency analysis will be performed.") + else: + logger.info("Statistial inefficiency analysis won't be performed.") G.analyze_alchemlyb(**kwargs) try: From c8a4ea43b54a0ea438463d7a4d5919e4ac02cf93 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Fri, 4 Jun 2021 04:58:08 -0700 Subject: [PATCH 13/20] print stride start and stop values within analysis functions --- mdpow/fep.py | 11 +++++++---- scripts/mdpow-solvationenergy | 29 +++++++++++++---------------- 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/mdpow/fep.py b/mdpow/fep.py index 5c5d3131..2d499475 100644 --- a/mdpow/fep.py +++ b/mdpow/fep.py @@ -950,6 +950,7 @@ def analyze(self, force=False, stride=None, autosave=True, ncorrel=25000): .. _p526: http://books.google.co.uk/books?id=XmyO2oRUg0cC&pg=PA526 """ stride = stride or self.stride + logger.info("Analysis stride is %s.",stride) if force or not self.has_dVdl(): try: @@ -1021,7 +1022,7 @@ def collect_alchemlyb(self, SI=True, start=0, stop=None, stride=None, autosave=T xvg_file = self.dgdl_xvg(self.wdir(component, l)) xvg_df = extract(xvg_file, T=self.Temperature).iloc[start:stop:stride] if SI: - logger.info("Performing statistical inefficiency analysis for window %s %s.", component, l) + logger.info("Performing statistical inefficiency analysis for window %s %04d" % (component, 1000 * l)) ts = _extract_dataframe(xvg_file).iloc[start:stop:stride] ts = pd.DataFrame({'time': ts.iloc[:,0], 'dhdl': ts.iloc[:,1]}) ts = ts.set_index('time') @@ -1044,6 +1045,10 @@ def analyze_alchemlyb(self, SI=True, start=0, stop=None, stride=None, force=Fals start = start or self.start stop = stop or self.stop + logger.info("Analysis stride is %s.",stride) + logger.info("Analysis starts from frame %s.",start) + logger.info("Analysis stops at frame %s.", stop) + if self.method in ['TI', 'BAR', 'MBAR']: estimator = self.estimators[self.method]['estimator'] else: @@ -1375,14 +1380,12 @@ def p_transfer(G1, G2, **kwargs): if kwargs['force'] or (not hasattr(G.results.DeltaA, 'Gibbs')): # write out the settings when the analysis is performed + logger.info("The solvent is %s .", G.solvent_type) logger.info("Estimator is %s.", estimator) logger.info("Free energy calculation method is %s.", G.method) - logger.info("Analysis stride is %s.", akw['stride']) if estimator == 'mdpow': G.analyze(**kwargs) elif estimator == 'alchemlyb': - logger.info("Analysis starts from frame %s.", G.start) - logger.info("Analysis stops at frame %s.", G.stop) if kwargs['SI']: logger.info("Statistial inefficiency analysis will be performed.") else: diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index 6b13e5cd..aa5f91d5 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -120,6 +120,9 @@ def run_gsolv(solvent, directory, **kwargs): logger.error(errmsg) raise ValueError(errmsg) + energyfile = kwargs.pop('energyfile', None) + plotfile = kwargs.pop('plotfile', None) + if not hasattr(gsolv, 'start'): gsolv.start = kwargs.pop('start', 0) if not hasattr(gsolv, 'stop'): @@ -129,27 +132,22 @@ def run_gsolv(solvent, directory, **kwargs): else: kwargs.setdefault('SI', gsolv.SI) - # do the analysis (only use keywords that are needed for FEP.analysis()) - akw = {'stride': kwargs.pop('stride', None), 'force': kwargs.pop('force', False)} - # make sure that we have data - if akw['force'] or (not hasattr(gsolv.results.DeltaA, 'Gibbs')): + if kwargs['force'] or (not hasattr(gsolv.results.DeltaA, 'Gibbs')): logger.info("[%(directory)s] Forcing (re)analysis of hydration free energy data", vars()) # write out the settings when the analysis is performed logger.info("Estimator is %s.", estimator) logger.info("Free energy calculation method is %s.", gsolv.method) - logger.info("Analysis stride is %s.", akw['stride']) if estimator == 'mdpow': + # do the analysis (only use keywords that are needed for FEP.analysis()) + akw = {'stride': kwargs.pop('stride', None), 'force': kwargs.pop('force', False)} gsolv.analyze(**akw) elif estimator == 'alchemlyb': - # the options only available for alchemlyb - logger.info("Analysis starts from frame %s.", gsolv.start) - logger.info("Analysis stops at frame %s.", gsolv.stop) if kwargs['SI']: logger.info("Statistial inefficiency analysis will be performed.") else: logger.info("Statistial inefficiency analysis won't be performed.") - gsolv.analyze_alchemlyb(**akw) + gsolv.analyze_alchemlyb(**kwargs) else: logger.info("[%(directory)s] Reading existing data from pickle file.", vars()) @@ -158,13 +156,12 @@ def run_gsolv(solvent, directory, **kwargs): if datstat(gsolv) == 'BAD': logger.warning("[%s] Possible file corruption: %s:%s", directory, solvent, datstat(gsolv)) - if kwargs.get('energyfile', None): + if energyfile: with mdpow.filelock.FileLock(kwargs['energyfile'], timeout=2): - with open(kwargs['energyfile'], mode='a') as out: + with open(energyfile, mode='a') as out: out.write(fmt_energy % (gsolv.summary(), directory)) logger.info("[%s] Wrote solvation energy terms to %r", directory, kwargs['energyfile']) - plotfile = kwargs.get('plotfile', None) if estimator == 'alchemlyb': plotfile = False if plotfile: @@ -286,10 +283,10 @@ if __name__ == "__main__": try: run_gsolv(opts.solvent, directory, plotfile=plotfile, - energyfile=realpath_or_none(opts.energyfile), - force=opts.force, stride=opts.stride, start=opts.start, - stop=opts.stop, estimator=opts.estimator, - method=opts.method, permissive=opts.permissive, SI=opts.SI) + energyfile=realpath_or_none(opts.energyfile), + force=opts.force, stride=opts.stride, start=opts.start, + stop=opts.stop, estimator=opts.estimator, + method=opts.method, permissive=opts.permissive, SI=opts.SI) except (OSError, IOError), err: logger.error("Running analysis in directory %r failed: %s", directory, str(err)) except Exception, err: From 5223563890565a34847776da9754b2ac45be7f4f Mon Sep 17 00:00:00 2001 From: VOD555 Date: Tue, 8 Jun 2021 02:35:01 -0700 Subject: [PATCH 14/20] use copy of kwargs in fep --- mdpow/fep.py | 26 +++++++++----------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/mdpow/fep.py b/mdpow/fep.py index 2d499475..492607cd 100644 --- a/mdpow/fep.py +++ b/mdpow/fep.py @@ -1360,18 +1360,19 @@ def p_transfer(G1, G2, **kwargs): logger.info("[%s] transfer free energy %s --> %s calculation", G1.molecule, G1.solvent_type, G2.solvent_type) for G in (G1, G2): + G_kwargs = kwargs.copy() # for fep files generated with old code which doesn't have these attributes if not hasattr(G, 'start'): - G.start = kwargs.pop('start', 0) + G.start = G_kwargs.pop('start', 0) if not hasattr(G, 'stop'): - G.stop = kwargs.pop('stop', None) + G.stop = G_kwargs.pop('stop', None) if not hasattr(G, 'SI'): - kwargs.setdefault('SI', True) + G_kwargs.setdefault('SI', True) else: - kwargs.setdefault('SI', G.SI) + G_kwargs.setdefault('SI', G.SI) # for this version. use the method given instead of the one in the input cfg file - G.method = kwargs.pop('method', 'MBAR') + G.method = G_kwargs.pop('method', 'MBAR') if estimator == 'mdpow': if G.method != "TI": errmsg = "Method %s is not implemented in MDPOW, use estimator='alchemlyb'" % G.method @@ -1384,23 +1385,14 @@ def p_transfer(G1, G2, **kwargs): logger.info("Estimator is %s.", estimator) logger.info("Free energy calculation method is %s.", G.method) if estimator == 'mdpow': - G.analyze(**kwargs) + G.analyze(**G_kwargs) elif estimator == 'alchemlyb': - if kwargs['SI']: + if G_kwargs['SI']: logger.info("Statistial inefficiency analysis will be performed.") else: logger.info("Statistial inefficiency analysis won't be performed.") - G.analyze_alchemlyb(**kwargs) + G.analyze_alchemlyb(**G_kwargs) - try: - G.results.DeltaA.Gibbs - G.logger_DeltaA0() - except (KeyError, AttributeError): # KeyError because results is a AttributeDict - logger.warn("Must analyze simulation because no hydration free energy data available...") - if estimator == 'mdpow': - G.analyze(**kwargs) - elif estimator == 'alchemlyb': - G.analyze_alchemlyb(**kwargs) # x.Gibbs are QuantityWithError so they do error propagation transferFE = G2.results.DeltaA.Gibbs - G1.results.DeltaA.Gibbs # note minus sign, with our convention of the free energy difference to be From 55e4ceb394ed68b9a6ca656bfc9ff66a989f57fc Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Tue, 8 Jun 2021 15:25:17 -0700 Subject: [PATCH 15/20] change SI option for mdpow-solvationenergy - changed `--noSI` to `--no-SI` - added `--SI` (the default) --- scripts/mdpow-pow | 2 +- scripts/mdpow-solvationenergy | 11 ++++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/scripts/mdpow-pow b/scripts/mdpow-pow index 2bbf01b2..3fd198b3 100755 --- a/scripts/mdpow-pow +++ b/scripts/mdpow-pow @@ -209,7 +209,7 @@ if __name__ == "__main__": parser.add_option('--noSI', dest='SI', action="store_false", help="Disable statistical inefficiency analysis" - "Statitical inefficiency analysis is performed by default when using" + "Statistical inefficiency analysis is performed by default when using" "alchemlyb estimators and is disabled when using mdpow estimator.") parser.add_option('-s', '--stride', dest="stride", type='int', help="Use every N-th datapoint from the original dV/dlambda data. " diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index aa5f91d5..9d8d2ac4 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -104,7 +104,7 @@ def run_gsolv(solvent, directory, **kwargs): return _datastatus[g.contains_corrupted_xvgs()] gsolv = load_gsolv(solvent, directory, permissive=kwargs.pop('permissive',False),) - logger.info("The solvent is %s .", solvent) + logger.info("The solvent is %s.", solvent) # for this version. use the method given instead of the one in the input cfg file estimator = kwargs.pop('estimator', 'alchemlyb') @@ -246,10 +246,15 @@ if __name__ == "__main__": parser.add_argument('--force', dest='force', default=False, action="store_true", help="force rereading all data ") - parser.add_argument('--noSI', dest='SI', + parser.add_argument('--SI', dest='SI', + action="store_true", default=True, + help="enable statistical inefficiency (SI) analysis. " + "Statistical inefficiency analysis is performed by default when using" + "alchemlyb estimators and is always disabled when using mdpow estimator.") + parser.add_argument('--no-SI', dest='SI', action="store_false", default=True, help="disable statistical inefficiency analysis. " - "Statitical inefficiency analysis is performed by default when using" + "Statistical inefficiency analysis is performed by default when using" "alchemlyb estimators and is disabled when using mdpow estimator. ") parser.add_argument('-s', '--stride', dest="stride", type=int, default=1, help="use every N-th datapoint from the original dV/dlambda data. ", From 2f4d54d875057823ed8a2cbd365bdd67c377edba Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Tue, 8 Jun 2021 15:27:53 -0700 Subject: [PATCH 16/20] Apply suggestions from code review --- mdpow/fep.py | 4 ++-- scripts/mdpow-solvationenergy | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/mdpow/fep.py b/mdpow/fep.py index 492607cd..753deb5a 100644 --- a/mdpow/fep.py +++ b/mdpow/fep.py @@ -1388,9 +1388,9 @@ def p_transfer(G1, G2, **kwargs): G.analyze(**G_kwargs) elif estimator == 'alchemlyb': if G_kwargs['SI']: - logger.info("Statistial inefficiency analysis will be performed.") + logger.info("Statistical inefficiency analysis will be performed.") else: - logger.info("Statistial inefficiency analysis won't be performed.") + logger.info("Statistical inefficiency analysis won't be performed.") G.analyze_alchemlyb(**G_kwargs) # x.Gibbs are QuantityWithError so they do error propagation diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index 9d8d2ac4..db275c09 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -144,9 +144,9 @@ def run_gsolv(solvent, directory, **kwargs): gsolv.analyze(**akw) elif estimator == 'alchemlyb': if kwargs['SI']: - logger.info("Statistial inefficiency analysis will be performed.") + logger.info("Statistical inefficiency analysis will be performed.") else: - logger.info("Statistial inefficiency analysis won't be performed.") + logger.info("Statistical inefficiency analysis won't be performed.") gsolv.analyze_alchemlyb(**kwargs) else: logger.info("[%(directory)s] Reading existing data from pickle file.", vars()) From 82527bb37d1c0f79a345c71052c95d86a36dbd96 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Tue, 8 Jun 2021 18:03:44 -0700 Subject: [PATCH 17/20] Raise error when SI=Ture and estimator=mdpow --- scripts/mdpow-solvationenergy | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index db275c09..5e9ed3ab 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -271,13 +271,14 @@ if __name__ == "__main__": "WRONG RESULTS CAN OCCUR! USE AT YOUR OWN RISK ") opts = parser.parse_args() - if len(opts.directory) == 0: - logger.fatal("A directory is required. See --help.") - sys.exit(1) - elif len(opts.directory) > 1 and not opts.plotfile.lower() in ('none', 'auto'): + if len(opts.directory) > 1 and not opts.plotfile.lower() in ('none', 'auto'): logger.fatal("Can only use --plotfile=None or --plotfile=auto with multiple directories.") sys.exit(1) + if opts.estimator == 'mdpow' and opt.SI: + logger.fatal("Statistical inefficiency analysis is only available for estimator 'alchemlyb'.") + sys.exit(1) + for directory in opts.directory: if not os.path.exists(directory): logger.warn("Directory %r not found, skipping...", directory) From b81ee0ac0bb5d4513e06a4f3adb4d617a11c2f2f Mon Sep 17 00:00:00 2001 From: VOD555 Date: Tue, 8 Jun 2021 18:23:00 -0700 Subject: [PATCH 18/20] fix error logger --- scripts/mdpow-solvationenergy | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index 5e9ed3ab..d4ce5cee 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -275,7 +275,7 @@ if __name__ == "__main__": logger.fatal("Can only use --plotfile=None or --plotfile=auto with multiple directories.") sys.exit(1) - if opts.estimator == 'mdpow' and opt.SI: + if opts.estimator == 'mdpow' and opts.SI: logger.fatal("Statistical inefficiency analysis is only available for estimator 'alchemlyb'.") sys.exit(1) From 650c3cf98d2e12a2cd3142e26b79399a986727e4 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Tue, 8 Jun 2021 18:31:58 -0700 Subject: [PATCH 19/20] pop keywords for final output before analysis --- scripts/mdpow-solvationenergy | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/mdpow-solvationenergy b/scripts/mdpow-solvationenergy index d4ce5cee..33bb334f 100755 --- a/scripts/mdpow-solvationenergy +++ b/scripts/mdpow-solvationenergy @@ -157,10 +157,10 @@ def run_gsolv(solvent, directory, **kwargs): logger.warning("[%s] Possible file corruption: %s:%s", directory, solvent, datstat(gsolv)) if energyfile: - with mdpow.filelock.FileLock(kwargs['energyfile'], timeout=2): + with mdpow.filelock.FileLock(energyfile, timeout=2): with open(energyfile, mode='a') as out: out.write(fmt_energy % (gsolv.summary(), directory)) - logger.info("[%s] Wrote solvation energy terms to %r", directory, kwargs['energyfile']) + logger.info("[%s] Wrote solvation energy terms to %r", directory, energyfile) if estimator == 'alchemlyb': plotfile = False From 202f4c97158ef47cdb6831dbcdb888544470bc14 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Tue, 8 Jun 2021 20:51:55 -0700 Subject: [PATCH 20/20] update CHANGES --- CHANGES | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CHANGES b/CHANGES index b24c60dd..311eeeb3 100644 --- a/CHANGES +++ b/CHANGES @@ -26,6 +26,14 @@ orbeckst, VOD555 invocation (#128) * Configuration setting mdrun.maxthreads now also applies to energy minimization with mdpow-equilibrium.py/mdpow.run() (#119) +* supported forcefield options in scripts (#123) +* supported multiple edr files input in fep.py (#132) +* supported alchemlyb estimators in analysis functions and scripts + (#133 #135) +* fixed: return fig in fep.py so plot function in scripts can access the + object (#129) +* support opls-aa tip4pd, amber tip4pew and cyclohexane, and charmm TIP5P + and cyclohexane solvent type (#141) 2017-05-02 0.6.1