From 7cdd5c8ef048676fea90435e46b947989627daa2 Mon Sep 17 00:00:00 2001 From: SimonBoothroyd Date: Tue, 16 Jan 2024 08:01:53 -0500 Subject: [PATCH] Add regression test analysis (#69) --- absolv/config.py | 26 - absolv/setup.py | 3 - absolv/tests/test_config.py | 28 +- docs/images/regression.svg | 4385 ++++++----------------------------- docs/reproducibility.md | 42 +- regression/README.md | 9 + regression/analyze.py | 111 + 7 files changed, 861 insertions(+), 3743 deletions(-) create mode 100644 regression/analyze.py diff --git a/absolv/config.py b/absolv/config.py index e4e730e..0d0af12 100644 --- a/absolv/config.py +++ b/absolv/config.py @@ -1,5 +1,4 @@ """Configure free energy calculations.""" -import math import typing import femto.md.config @@ -355,28 +354,3 @@ class Result(pydantic.BaseModel): dg_std_solvent_b: OpenMMQuantity[KCAL_MOL] = pydantic.Field( description="The standard error in ``dg_solvent_b``." ) - - @property - def dg(self) -> openmm.unit.Quantity: - """The change in free energy of transferring the solute from *solvent-a* to - *solvent-b* in units of kT.""" - return self.dg_solvent_b - self.dg_solvent_a - - @property - def dg_std(self) -> openmm.unit.Quantity: - """The standard error in ``ddg``.""" - - std = math.sqrt( - self.dg_std_solvent_a.value_in_unit(KCAL_MOL) ** 2 - + self.dg_std_solvent_b.value_in_unit(KCAL_MOL) ** 2 - ) - return std * KCAL_MOL - - def __str__(self): - return ( - f"ΔG a->b={self.dg.value_in_unit(KCAL_MOL):.3f} kcal/mol " - f"ΔG a->b std={self.dg_std.value_in_unit(KCAL_MOL):.3f} kcal/mol" - ) - - def __repr__(self): - return f"{self.__repr_name__()}({self.__str__()})" diff --git a/absolv/setup.py b/absolv/setup.py index 162e880..7b7baf7 100644 --- a/absolv/setup.py +++ b/absolv/setup.py @@ -118,9 +118,6 @@ def setup_system( tolerance: The minimum spacing between molecules during packing in units compatible with angstroms. - Raises: - * PACKMOLRuntimeError - Returns: A topology containing the molecules the coordinates were generated for and a unit [A] wrapped numpy array of coordinates with shape=(n_atoms, 3). diff --git a/absolv/tests/test_config.py b/absolv/tests/test_config.py index 3fbd2cb..64e8851 100644 --- a/absolv/tests/test_config.py +++ b/absolv/tests/test_config.py @@ -1,10 +1,7 @@ -import math - -import openmm.unit import pydantic import pytest -from absolv.config import EquilibriumProtocol, Result, System +from absolv.config import EquilibriumProtocol, System class TestSystem: @@ -80,26 +77,3 @@ def test_validate_lambda_lengths(self, lambda_sterics, lambda_electrostatics): lambda_sterics=lambda_sterics, lambda_electrostatics=lambda_electrostatics, ) - - -class TestResult: - @pytest.fixture - def mock_result(self): - return Result( - dg_solvent_a=1.0 * openmm.unit.kilocalorie_per_mole, - dg_std_solvent_a=0.2 * openmm.unit.kilocalorie_per_mole, - dg_solvent_b=3.0 * openmm.unit.kilocalorie_per_mole, - dg_std_solvent_b=0.4 * openmm.unit.kilocalorie_per_mole, - ) - - def test_dg(self, mock_result): - assert mock_result.dg == pytest.approx(2.0) * openmm.unit.kilocalorie_per_mole - assert ( - mock_result.dg_std - == pytest.approx(math.sqrt(0.2**2 + 0.4**2)) - * openmm.unit.kilocalorie_per_mole - ) - - def test_repr(self, mock_result): - assert repr(mock_result) is not None - assert str(mock_result) is not None diff --git a/docs/images/regression.svg b/docs/images/regression.svg index 3f439f8..fe34bf2 100644 --- a/docs/images/regression.svg +++ b/docs/images/regression.svg @@ -1,186 +1,152 @@ - + - + - 2021-12-23T15:25:27.606607 + 2024-01-16T07:38:02.650012 image/svg+xml - Matplotlib v3.4.3, https://matplotlib.org/ + Matplotlib v3.8.2, https://matplotlib.org/ - + - +" style="fill: #ffffff"/> - +" style="fill: #ffffff"/> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + - + - +" style="stroke: #000000; stroke-width: 0.8"/> - + - - - - - + + + + + + + + + - + - - - - - - - - - - - + - - - - - + + + + + + + + - - - - - - - + + - + - - - + + + - - - - + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +" transform="scale(0.015625)"/> - - - - - - - - - - - - - - - - - - - - + - - + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + - - - - - - + + + - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + - + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +" transform="scale(0.015625)"/> + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - + + + + + + - + - - - + + + - - - + - - + + - + - - - + + + - - - + - - + + - + - - - + + + - - - + - - + + - + - - - - - - - + + + + - - + + - + - - - - - - - + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + - - - + - + - + - + - +" transform="scale(0.015625)"/> - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - +" transform="scale(0.015625)"/> - - - - - - - + - - + + - - - - - - - - - - - - - - - + + - - + + + + + - + - + - - - - - - - - - - - - - - - - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - + + diff --git a/docs/reproducibility.md b/docs/reproducibility.md index 9eff500..caeaf1b 100644 --- a/docs/reproducibility.md +++ b/docs/reproducibility.md @@ -1,29 +1,31 @@ # Reproducibility -While every effort has been made to ensure the 'correctness and reproducibility' of any results computed using this -framework, achieving consistent free energies between different frameworks and simulation engines has been notoriously -tricky [@loeffler2018reproducibility]. +While every effort has been made to ensure the 'correctness and reproducibility' of any +results computed using this framework, achieving consistent free energies between +different frameworks and simulation engines has been notoriously tricky +[@loeffler2018reproducibility]. -In an attempt to ensure that this framework remains at least self-consistent between versions, and as consistent as -possible with other packages, a suite of regression tests are provided in the main [GitHub repository](https://github.com/SimonBoothroyd/absolv/). +In an attempt to ensure that this framework remains at least self-consistent between +versions, and as consistent as possible with other packages, a suite of regression tests +are provided in the main [GitHub repository](https://github.com/SimonBoothroyd/absolv/). -These include tests to ensure that: - -- the energies of systems that been alchemically modified are consistent with the original system -- systems that contain custom non-bonded force are alchemically transformed in the exact same was as systems that - contain the built-in OpenMM LJ non-bonded force - -and most importantly, that computing the free energies using both the 'equilibrium' and 'non-equilibrium' methods -supported in this framework are in agreement amongst themselves, with values computed using [Yank](https://github.com/choderalab/yank), and with the -GROMACS values reported by Loeffler *et al* [@loeffler2018reproducibility]. +These include tests to ensure that computing the free energies using both the +'equilibrium' and 'non-equilibrium' methods supported in this framework are in agreement +amongst themselves, and with the GROMACS values reported by Loeffler *et al* +[@loeffler2018reproducibility]. ## Regression Results -The results of running the free energy comparisons using the latest version of the framework are shown below: - -![](images/regression.svg){.align-center width="80.0%"} - -Here the error bars show the standard deviation computed across three -replicas. +The results of running the free energy comparisons using the latest version of the +framework are shown below: + +
+ ![Regression Results](images/regression.svg){ width="100.0%" } +
The average hydration free energies of methane, methanol, toluene, and + 2-cyclopentanylindole computed using the 'equilibrium' and + 'non-equilibrium' methods across three replicas. The error bars show the + standard deviation computed across three replicas. +
+
\bibliography diff --git a/regression/README.md b/regression/README.md index 9bd5889..5a20733 100644 --- a/regression/README.md +++ b/regression/README.md @@ -29,3 +29,12 @@ python run.py --solute methanol \ ``` The results will be written to `results/{timestamp}/{method}-{solute}-{replica}/results.json`. + +The results can be plotted using: + +```shell +python analyze.py --result-dir results/2024-01-15--12-00-40 \ + --result-dir results/2024-01-15--18-20-15 \ + --result-dir ... \ + --output results/1.0.0.svg +``` diff --git a/regression/analyze.py b/regression/analyze.py new file mode 100644 index 0000000..9b654a6 --- /dev/null +++ b/regression/analyze.py @@ -0,0 +1,111 @@ +import collections +import json +import logging +import pathlib + +import click +import numpy +import openmm.unit +from matplotlib import pyplot +from run import DEFAULT_METHODS, DEFAULT_SOLUTES + +import absolv.config + +KCAL = openmm.unit.kilocalorie_per_mole + + +def _compute_dg(result: absolv.config.Result) -> float: + # for the regression test we set up solvent A as vacuum and solvent B as + # explicit solvent. + # in both 'solvents' the alchemical path goes from molecule switched on to + # molecule switched off. + return (result.dg_solvent_a + -result.dg_solvent_b).value_in_unit(KCAL) + + +def _parse_results( + result_dirs: list[pathlib.Path], +) -> dict[str, dict[str, tuple[float, float]]]: + results = collections.defaultdict(dict) + + for method in DEFAULT_METHODS: + for solute in DEFAULT_SOLUTES: + replica_paths = [ + result_path + for result_dir in result_dirs + for result_path in result_dir.glob(f"{method}-{solute}-*/result.json") + ] + + if len(replica_paths) == 0: + continue + + replicas = [ + absolv.config.Result.model_validate_json(result_path.read_text()) + for result_path in replica_paths + ] + + dg = numpy.array([_compute_dg(replica) for replica in replicas]) + results[method][solute] = numpy.mean(dg), numpy.std(dg) + + return {**results} + + +@click.command() +@click.option( + "--result-dir", + "result_dirs", + multiple=True, + type=click.Path( + exists=True, file_okay=False, dir_okay=True, path_type=pathlib.Path + ), +) +@click.option( + "--output", + "output_path", + type=click.Path( + exists=False, file_okay=True, dir_okay=False, path_type=pathlib.Path + ), + required=True, +) +def main(result_dirs: list[pathlib.Path], output_path: pathlib.Path): + logging.basicConfig(level=logging.INFO) + + reference = json.loads(pathlib.Path("results/reference.json").read_text()) + results = _parse_results(result_dirs) + + figure, axis = pyplot.subplots(1, 1, figsize=(4, 4)) + + for method in results: + solutes = sorted(results[method]) + + x = [reference[solute]["value"] for solute in solutes] + x_std = [reference[solute]["std_error"] for solute in solutes] + + y = [results[method][solute][0] for solute in solutes] + y_std = [results[method][solute][1] for solute in solutes] + + axis.errorbar( + x, y, xerr=x_std, yerr=y_std, label=method, marker="x", linestyle="none" + ) + + lims = [ + numpy.min([axis.get_xlim(), axis.get_ylim()]), + numpy.max([axis.get_xlim(), axis.get_ylim()]), + ] + + # now plot both limits against eachother + axis.plot(lims, lims, "k-", alpha=0.25, zorder=0) + axis.set_aspect("equal") + axis.set_xlim(lims) + axis.set_ylim(lims) + + pyplot.xlabel(r"Reference $\Delta G$ [kcal/mol]") + pyplot.ylabel(r"Estimated $\Delta G$ [kcal/mol]") + + pyplot.legend() + + pyplot.tight_layout() + pyplot.savefig(output_path) + + +if __name__ == "__main__": + main()