From 97a0a289cfe3c76041704a0ecd29e19dfb75f60f Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Thu, 22 Oct 2026 11:40:46 +0100 Subject: [PATCH 01/11] exploring rbfe settings --- cookbook/exploring_rbfe_settings.ipynb | 292 +++++++++++++++++++++++++ 1 file changed, 292 insertions(+) create mode 100644 cookbook/exploring_rbfe_settings.ipynb diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb new file mode 100644 index 0000000..d4d6b94 --- /dev/null +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -0,0 +1,292 @@ +{ + "cells": [ + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "# Exploring HybridTopologyProtocol settings\n", + "\n", + "The `HybridTopologyProtocol` exposes a number of simulation settings that can be used to control the behavior of an alchemical transformation. In this notebook, we will explore changes to common settings, explaining when you might want to\n", + "consider changing them, the expected effect of changing them, and how to change them in practice.\n", + "\n", + "## Inspecting the default settings\n", + "\n", + "The default settings for the `HybridTopologyProtocol` are designed to be a good starting point for most alchemical transformations. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision.\n", + "\n", + "The default settings can be inspected by:" + ], + "id": "a3f965b58be21979" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-21T17:14:24.838663Z", + "start_time": "2026-10-21T17:14:24.718322Z" + } + }, + "cell_type": "code", + "source": [ + "from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol\n", + "\n", + "settings = RelativeHybridTopologyProtocol.default_settings()" + ], + "id": "5690f07a5cd27245", + "outputs": [], + "execution_count": 10 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "For a complete explanation of each setting and how to customize the object see [Choose and Configure a Protocol](https://docs.openfree.energy/en/latest/cookbook/choose_protocol.html).\n", + "\n", + "The main parameters that determine the amount of sampling performed and the associated cost for the hybrid topology protocol are:\n", + "\n", + "- **Lambda windows**: the number of discrete alchemical states simulated during the transformation.\n", + "- **Production length**: the length of the production simulation performed at each lambda window.\n", + "- **Equilibration length**: the length of the equilibration simulation performed at each lambda window.\n", + "- **Protocol repeats**: the number of independent executions of the protocol.\n", + "\n", + "These parameters can be adjusted to increase the amount of sampling performed, which can improve convergence and precision, but will also increase the computational cost of the simulation. The following example retrieves them from the default settings:" + ], + "id": "94c914fc75cd8fcb" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-21T17:14:25.858033Z", + "start_time": "2026-10-21T17:14:25.812345Z" + } + }, + "cell_type": "code", + "source": [ + "n_repeats = settings.protocol_repeats\n", + "n_windows = settings.lambda_settings.lambda_windows\n", + "production_length = settings.simulation_settings.production_length\n", + "equilibration_length = settings.simulation_settings.equilibration_length\n", + "\n", + "print(f\"Lambda windows: {n_windows}\")\n", + "print(f\"Production length per state: {production_length}\")\n", + "print(f\"Equilibration length per state: {equilibration_length}\")\n", + "print(f\"Independent repeats: {n_repeats}\")" + ], + "id": "3c2f7354aadac1e8", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Lambda windows: 11\n", + "Production length per state: 5.0 nanosecond\n", + "Equilibration length per state: 1.0 nanosecond\n", + "Independent repeats: 3\n" + ] + } + ], + "execution_count": 11 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## Estimating the computational cost of a protocol\n", + "\n", + "A useful first approximation of the computational cost of a protocol execution is the total simulation time in nanoseconds across all lambda windows and repeats. This can be calculated as:\n", + "\n", + "$$T_{total} = N_{repeats} \\times N_{\\lambda windows} \\times (T_{production} + T_{equilibration})$$\n", + "\n", + "So for the default settings the total simulation time is:" + ], + "id": "704d5dc6e6f895f5" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-21T17:14:27.457661Z", + "start_time": "2026-10-21T17:14:27.409115Z" + } + }, + "cell_type": "code", + "source": [ + "total_sim_time = n_repeats * n_windows * (production_length + equilibration_length)\n", + "print(f\"Total simulation time: {total_sim_time} ns\")" + ], + "id": "fd81e13e2c34684d", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total simulation time: 198.0 nanosecond ns\n" + ] + } + ], + "execution_count": 12 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "The HybridTopologyProtocol simulates the complex and solvent legs of a transformation independently, so if both use the same settings the total cost of the protocol is double the above estimate. It is useful to consider this estimate when deciding how to adjust the settings for a given transformation.\n", + "\n", + "# Increasing simulation length\n", + "\n", + "This section describes how to increase the amount of production simulation performed at each lambda window and estimate the increase in computational cost.\n", + "\n", + "## When to increase simulation length\n", + "\n", + "Consider increasing the simulation length when:\n", + "\n", + "- Free-energy estimates differ substantially between independent repeats of the protocol.\n", + "- Forward and reverse free-energy estimates disagree by more than 1kbT 50% of the simulation length and their uncertainties do not overlap.\n", + "- A large number of samples are discarded during analysis due to the automated equilibration and decorrelation methods.\n", + "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows.\n", + "- The transformation involves large structural changes due to for example changes in formal charge.\n", + "\n", + "This information is available through the automated analysis performed by the `HybridTopologyProtocol` see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for more information, it is recommended to review the analysis metrics before deciding to increase the simulation length.\n", + "\n", + "## When not to increase simulation length\n", + "\n", + "Increasing the simulation length is not always the best solution. For example poor overlap between neighboring lambda windows may require changes to the number of lambda windows or the spacing between them.\n", + "\n", + "## Changing the production length\n", + "\n", + "The following example shows how to change the production length for a `HybridTopologyProtocol` object. In this example, we will double the production length from `5 ns` to `10 ns` per lambda window. Note the use of units from the `openff.units` package to specify the value is in nanoseconds. We will then calculate the new total simulation time and compare it to the default settings.\n" + ], + "id": "3ab10c2ecebbd3df" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-22T10:36:48.805701Z", + "start_time": "2026-10-22T10:36:48.498120Z" + } + }, + "cell_type": "code", + "source": [ + "from openff.units import unit\n", + "\n", + "long_settings = RelativeHybridTopologyProtocol.default_settings()\n", + "\n", + "# Increase the production length to 10 ns\n", + "long_settings.simulation_settings.production_length = 10 * unit.nanoseconds\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "protocol = RelativeHybridTopologyProtocol(settings=long_settings)\n", + "\n", + "print(f\"Updated production length: {protocol.settings.simulation_settings.production_length}\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_total_sim_time = n_repeats * n_windows * (protocol.settings.simulation_settings.production_length + protocol.settings.simulation_settings.equilibration_length)\n", + "print(f\"New total simulation time: {new_total_sim_time} ns\")\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_total_sim_time - total_sim_time) / total_sim_time) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "2cf9b3c7030df36f", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated production length: 10 nanosecond\n", + "New total simulation time: 363.0 nanosecond ns\n", + "Percentage increase in total simulation time: 83.33%\n" + ] + } + ], + "execution_count": 22 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "# Increasing the number of lambda windows\n", + "\n", + "Increasing the number of lambda windows introduces additional intermediate alchemical states along the transformation which can make the pathway easier to traverse and improve the overlap between neighboring states. This section describes how to increase the number of lambda windows and estimate the increase in computational cost.\n", + "\n", + "## When to increase the number of lambda windows\n", + "\n", + "Consider increasing the number of lambda windows **before** running the simulation when:\n", + "\n", + "- The transformation involves a change in net formal charge or substantial redistribution of charge.\n", + "\n", + "Consider increasing the number of lambda windows **after** running the simulation when:\n", + "\n", + "- Neighboring lambda windows have poor phase-space overlap or replica mixing statistics.\n", + "\n", + "## When not to increase the number of lambda windows\n", + "\n", + "Increasing the number of lambda windows is not always the best solution, slow changes in ligand conformation or binding site rearrangements may require longer simulation lengths to achieve convergence. In these cases, increasing the number of lambda windows may not improve convergence and will increase the computational cost of the protocol.\n", + "\n", + "## Changing the number of lambda windows\n", + "\n", + "The following example shows how to change the number of lambda windows for a `HybridTopologyProtocol` object. In this example, we will double the number of lambda windows from `11` to `22`. We will then calculate the new total simulation time and compare it to the default settings. Note the `lambda_windows` setting only accepts integer values and generates a uniform spacing of lambda values between `0` and `1`.\n" + ], + "id": "3d6bb208a91ce6d0" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-22T10:38:21.632585Z", + "start_time": "2026-10-22T10:38:21.575761Z" + } + }, + "cell_type": "code", + "source": [ + "more_window_settings = RelativeHybridTopologyProtocol.default_settings()\n", + "\n", + "# Increase the number of lambda windows to 22\n", + "more_window_settings.lambda_settings.lambda_windows = 22\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "protocol = RelativeHybridTopologyProtocol(settings=more_window_settings)\n", + "\n", + "print(f\"Updated number of lambda windows: {protocol.settings.lambda_settings.lambda_windows}\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_total_sim_time = n_repeats * protocol.settings.lambda_settings.lambda_windows * (protocol.settings.simulation_settings.production_length + protocol.settings.simulation_settings.equilibration_length)\n", + "print(f\"New total simulation time: {new_total_sim_time} ns\")\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_total_sim_time - total_sim_time) / total_sim_time) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "a8a0a2b5ba644eb4", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated number of lambda windows: 22\n", + "New total simulation time: 396.0 nanosecond ns\n", + "Percentage increase in total simulation time: 100.00%\n" + ] + } + ], + "execution_count": 23 + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 5c447f3fb8063aebfdb9c33e2f56aca25e5a0f3d Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Thu, 22 Oct 2026 11:46:42 +0100 Subject: [PATCH 02/11] small format update --- cookbook/exploring_rbfe_settings.ipynb | 40 ++++++++++++++------------ 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index d4d6b94..241985d 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -6,6 +6,7 @@ "source": [ "# Exploring HybridTopologyProtocol settings\n", "\n", + "\n", "The `HybridTopologyProtocol` exposes a number of simulation settings that can be used to control the behavior of an alchemical transformation. In this notebook, we will explore changes to common settings, explaining when you might want to\n", "consider changing them, the expected effect of changing them, and how to change them in practice.\n", "\n", @@ -20,8 +21,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-21T17:14:24.838663Z", - "start_time": "2026-10-21T17:14:24.718322Z" + "end_time": "2026-10-22T10:46:22.362791Z", + "start_time": "2026-10-22T10:46:22.275865Z" } }, "cell_type": "code", @@ -32,7 +33,7 @@ ], "id": "5690f07a5cd27245", "outputs": [], - "execution_count": 10 + "execution_count": 24 }, { "metadata": {}, @@ -54,8 +55,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-21T17:14:25.858033Z", - "start_time": "2026-10-21T17:14:25.812345Z" + "end_time": "2026-10-22T10:46:22.931214Z", + "start_time": "2026-10-22T10:46:22.867745Z" } }, "cell_type": "code", @@ -83,7 +84,7 @@ ] } ], - "execution_count": 11 + "execution_count": 25 }, { "metadata": {}, @@ -102,8 +103,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-21T17:14:27.457661Z", - "start_time": "2026-10-21T17:14:27.409115Z" + "end_time": "2026-10-22T10:46:23.520808Z", + "start_time": "2026-10-22T10:46:23.480068Z" } }, "cell_type": "code", @@ -121,7 +122,7 @@ ] } ], - "execution_count": 12 + "execution_count": 26 }, { "metadata": {}, @@ -135,13 +136,16 @@ "\n", "## When to increase simulation length\n", "\n", - "Consider increasing the simulation length when:\n", + "Consider increasing the simulation length **before** running the simulation when:\n", + "\n", + "- The transformation involves large structural changes due to for example changes in formal charge, you may also need to consider increasing the number of lambda windows see the next section for more information.\n", + "\n", + "Consider increasing the simulation length **after** running the simulation when:\n", "\n", "- Free-energy estimates differ substantially between independent repeats of the protocol.\n", "- Forward and reverse free-energy estimates disagree by more than 1kbT 50% of the simulation length and their uncertainties do not overlap.\n", "- A large number of samples are discarded during analysis due to the automated equilibration and decorrelation methods.\n", "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows.\n", - "- The transformation involves large structural changes due to for example changes in formal charge.\n", "\n", "This information is available through the automated analysis performed by the `HybridTopologyProtocol` see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for more information, it is recommended to review the analysis metrics before deciding to increase the simulation length.\n", "\n", @@ -158,8 +162,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T10:36:48.805701Z", - "start_time": "2026-10-22T10:36:48.498120Z" + "end_time": "2026-10-22T10:46:24.174430Z", + "start_time": "2026-10-22T10:46:24.112991Z" } }, "cell_type": "code", @@ -196,7 +200,7 @@ ] } ], - "execution_count": 22 + "execution_count": 27 }, { "metadata": {}, @@ -210,7 +214,7 @@ "\n", "Consider increasing the number of lambda windows **before** running the simulation when:\n", "\n", - "- The transformation involves a change in net formal charge or substantial redistribution of charge.\n", + "- The transformation involves a change in net formal charge or substantial redistribution of charge, you may also need to consider increasing the simulation length see the previous section for more information.\n", "\n", "Consider increasing the number of lambda windows **after** running the simulation when:\n", "\n", @@ -229,8 +233,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T10:38:21.632585Z", - "start_time": "2026-10-22T10:38:21.575761Z" + "end_time": "2026-10-22T10:46:24.818876Z", + "start_time": "2026-10-22T10:46:24.776710Z" } }, "cell_type": "code", @@ -265,7 +269,7 @@ ] } ], - "execution_count": 23 + "execution_count": 28 } ], "metadata": { From 4653a6e4f547f4c41f089af4d383babaf5c180c8 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Thu, 22 Oct 2026 12:33:06 +0100 Subject: [PATCH 03/11] clean up units and some text --- cookbook/exploring_rbfe_settings.ipynb | 49 +++++++++++++------------- 1 file changed, 24 insertions(+), 25 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index 241985d..01076f2 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -7,12 +7,11 @@ "# Exploring HybridTopologyProtocol settings\n", "\n", "\n", - "The `HybridTopologyProtocol` exposes a number of simulation settings that can be used to control the behavior of an alchemical transformation. In this notebook, we will explore changes to common settings, explaining when you might want to\n", - "consider changing them, the expected effect of changing them, and how to change them in practice.\n", + "The `HybridTopologyProtocol` exposes a number of simulation settings that can effect the convergence properties, accuracy and precision of an alchemical transformation. In this notebook, we will explore changes to common settings, explaining when you might want to consider changing them, the expected effect of changing them, and how to change them in practice.\n", "\n", "## Inspecting the default settings\n", "\n", - "The default settings for the `HybridTopologyProtocol` are designed to be a good starting point for most alchemical transformations. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision.\n", + "The default settings for the `HybridTopologyProtocol` are designed to be a good starting point for most alchemical transformations. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", "\n", "The default settings can be inspected by:" ], @@ -48,15 +47,15 @@ "- **Equilibration length**: the length of the equilibration simulation performed at each lambda window.\n", "- **Protocol repeats**: the number of independent executions of the protocol.\n", "\n", - "These parameters can be adjusted to increase the amount of sampling performed, which can improve convergence and precision, but will also increase the computational cost of the simulation. The following example retrieves them from the default settings:" + "The appropriate changes to these settings are transformation dependent with no single setting guarantying sufficient convergence or accuracy in general. These settings are good starting point to investigate the effect of additional sampling on the transformation under study however, but will also increase the computational cost of the simulation. The following example retrieves them from the default settings:" ], "id": "94c914fc75cd8fcb" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T10:46:22.931214Z", - "start_time": "2026-10-22T10:46:22.867745Z" + "end_time": "2026-10-22T11:21:12.003187Z", + "start_time": "2026-10-22T11:21:11.965473Z" } }, "cell_type": "code", @@ -67,8 +66,8 @@ "equilibration_length = settings.simulation_settings.equilibration_length\n", "\n", "print(f\"Lambda windows: {n_windows}\")\n", - "print(f\"Production length per state: {production_length}\")\n", - "print(f\"Equilibration length per state: {equilibration_length}\")\n", + "print(f\"Production length per state: {production_length.m} ns\")\n", + "print(f\"Equilibration length per state: {equilibration_length.m} ns\")\n", "print(f\"Independent repeats: {n_repeats}\")" ], "id": "3c2f7354aadac1e8", @@ -78,13 +77,13 @@ "output_type": "stream", "text": [ "Lambda windows: 11\n", - "Production length per state: 5.0 nanosecond\n", - "Equilibration length per state: 1.0 nanosecond\n", + "Production length per state: 5.0 ns\n", + "Equilibration length per state: 1.0 ns\n", "Independent repeats: 3\n" ] } ], - "execution_count": 25 + "execution_count": 32 }, { "metadata": {}, @@ -103,14 +102,14 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T10:46:23.520808Z", - "start_time": "2026-10-22T10:46:23.480068Z" + "end_time": "2026-10-22T11:20:41.651615Z", + "start_time": "2026-10-22T11:20:41.533380Z" } }, "cell_type": "code", "source": [ "total_sim_time = n_repeats * n_windows * (production_length + equilibration_length)\n", - "print(f\"Total simulation time: {total_sim_time} ns\")" + "print(f\"Total simulation time: {total_sim_time.m} ns\")" ], "id": "fd81e13e2c34684d", "outputs": [ @@ -118,11 +117,11 @@ "name": "stdout", "output_type": "stream", "text": [ - "Total simulation time: 198.0 nanosecond ns\n" + "Total simulation time: 198.0 ns\n" ] } ], - "execution_count": 26 + "execution_count": 29 }, { "metadata": {}, @@ -143,9 +142,9 @@ "Consider increasing the simulation length **after** running the simulation when:\n", "\n", "- Free-energy estimates differ substantially between independent repeats of the protocol.\n", - "- Forward and reverse free-energy estimates disagree by more than 1kbT 50% of the simulation length and their uncertainties do not overlap.\n", + "- Forward and reverse free-energy estimates disagree by more than 1kbT at 50% of the simulation length and their uncertainties do not overlap.\n", "- A large number of samples are discarded during analysis due to the automated equilibration and decorrelation methods.\n", - "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows.\n", + "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows, you may also need to consider increasing the number of lambda windows see the next section for more information.\n", "\n", "This information is available through the automated analysis performed by the `HybridTopologyProtocol` see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for more information, it is recommended to review the analysis metrics before deciding to increase the simulation length.\n", "\n", @@ -162,8 +161,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T10:46:24.174430Z", - "start_time": "2026-10-22T10:46:24.112991Z" + "end_time": "2026-10-22T11:21:00.199353Z", + "start_time": "2026-10-22T11:21:00.150806Z" } }, "cell_type": "code", @@ -178,11 +177,11 @@ "# Create a new protocol object with the updated settings\n", "protocol = RelativeHybridTopologyProtocol(settings=long_settings)\n", "\n", - "print(f\"Updated production length: {protocol.settings.simulation_settings.production_length}\")\n", + "print(f\"Updated production length: {protocol.settings.simulation_settings.production_length.m} ns\")\n", "\n", "# Calculate the new total simulation time\n", "new_total_sim_time = n_repeats * n_windows * (protocol.settings.simulation_settings.production_length + protocol.settings.simulation_settings.equilibration_length)\n", - "print(f\"New total simulation time: {new_total_sim_time} ns\")\n", + "print(f\"New total simulation time: {new_total_sim_time.m} ns\")\n", "\n", "# Print the percentage increase in total simulation time\n", "percentage_increase = ((new_total_sim_time - total_sim_time) / total_sim_time) * 100\n", @@ -194,13 +193,13 @@ "name": "stdout", "output_type": "stream", "text": [ - "Updated production length: 10 nanosecond\n", - "New total simulation time: 363.0 nanosecond ns\n", + "Updated production length: 10 ns\n", + "New total simulation time: 363.0 ns\n", "Percentage increase in total simulation time: 83.33%\n" ] } ], - "execution_count": 27 + "execution_count": 31 }, { "metadata": {}, From ded8d29c3bc6bfa46ef886e9ccc3b8f4758bbd19 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Thu, 22 Oct 2026 12:35:18 +0100 Subject: [PATCH 04/11] more intro --- cookbook/exploring_rbfe_settings.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index 01076f2..8abaee5 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -11,7 +11,7 @@ "\n", "## Inspecting the default settings\n", "\n", - "The default settings for the `HybridTopologyProtocol` are designed to be a good starting point for most alchemical transformations. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", + "The default settings for the `HybridTopologyProtocol` are designed to be a good starting point for most alchemical transformations with the goal of obtaining reliable estimates with sufficient precision for decision-making at a reasonable computational cost. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", "\n", "The default settings can be inspected by:" ], From 2987975358090c27f32cff3ef87feecc25feaf57 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Thu, 22 Oct 2026 14:21:25 +0100 Subject: [PATCH 05/11] more clean up --- cookbook/exploring_rbfe_settings.ipynb | 42 ++++++++++++++++---------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index 8abaee5..767d5b4 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -7,11 +7,11 @@ "# Exploring HybridTopologyProtocol settings\n", "\n", "\n", - "The `HybridTopologyProtocol` exposes a number of simulation settings that can effect the convergence properties, accuracy and precision of an alchemical transformation. In this notebook, we will explore changes to common settings, explaining when you might want to consider changing them, the expected effect of changing them, and how to change them in practice.\n", + "The `RelativeHybridTopologyProtocol` exposes a number of simulation settings that can affect the convergence properties, accuracy and precision of an alchemical transformation. In this notebook, we will explore changes to common settings, advising when you might want to consider changing them and how to change them in practice. For a detailed explanation of the protocol and its limitations see the [protocol documentation](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#relative-hybrid-topology-protocol) and references therein.\n", "\n", "## Inspecting the default settings\n", "\n", - "The default settings for the `HybridTopologyProtocol` are designed to be a good starting point for most alchemical transformations with the goal of obtaining reliable estimates with sufficient precision for decision-making at a reasonable computational cost. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", + "The default settings for the `RelativeHybridTopologyProtocol` are designed to be a good starting point for most alchemical transformations with the goal of obtaining reliable estimates with sufficient precision for decision-making at a reasonable computational cost. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", "\n", "The default settings can be inspected by:" ], @@ -47,15 +47,15 @@ "- **Equilibration length**: the length of the equilibration simulation performed at each lambda window.\n", "- **Protocol repeats**: the number of independent executions of the protocol.\n", "\n", - "The appropriate changes to these settings are transformation dependent with no single setting guarantying sufficient convergence or accuracy in general. These settings are good starting point to investigate the effect of additional sampling on the transformation under study however, but will also increase the computational cost of the simulation. The following example retrieves them from the default settings:" + "The appropriate changes to these settings are transformation dependent with no single setting guaranteeing sufficient convergence or accuracy in general. These settings provide a good starting point to investigate the affect of additional sampling on the transformation under study however, but will also increase the computational cost of the simulation. The following example retrieves them from the default settings:" ], "id": "94c914fc75cd8fcb" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T11:21:12.003187Z", - "start_time": "2026-10-22T11:21:11.965473Z" + "end_time": "2026-10-22T13:12:34.530752Z", + "start_time": "2026-10-22T13:12:34.485473Z" } }, "cell_type": "code", @@ -83,7 +83,7 @@ ] } ], - "execution_count": 32 + "execution_count": 35 }, { "metadata": {}, @@ -127,7 +127,7 @@ "metadata": {}, "cell_type": "markdown", "source": [ - "The HybridTopologyProtocol simulates the complex and solvent legs of a transformation independently, so if both use the same settings the total cost of the protocol is double the above estimate. It is useful to consider this estimate when deciding how to adjust the settings for a given transformation.\n", + "The `RelativeHybridTopologyProtocol` simulates the complex and solvent legs of a transformation independently, so if both use the same settings the total cost of the protocol is double the above estimate. It is useful to consider this estimate when deciding how to adjust the settings for a given transformation.\n", "\n", "# Increasing simulation length\n", "\n", @@ -142,11 +142,11 @@ "Consider increasing the simulation length **after** running the simulation when:\n", "\n", "- Free-energy estimates differ substantially between independent repeats of the protocol.\n", - "- Forward and reverse free-energy estimates disagree by more than 1kbT at 50% of the simulation length and their uncertainties do not overlap.\n", + "- Forward and reverse free-energy estimates disagree by more than 1kbT at over 50% of the simulation length and their uncertainties do not overlap.\n", "- A large number of samples are discarded during analysis due to the automated equilibration and decorrelation methods.\n", "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows, you may also need to consider increasing the number of lambda windows see the next section for more information.\n", "\n", - "This information is available through the automated analysis performed by the `HybridTopologyProtocol` see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for more information, it is recommended to review the analysis metrics before deciding to increase the simulation length.\n", + "This information is available through the automated analysis performed by the `RelativeHybridTopologyProtocol` see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for more information, it is recommended to review the analysis metrics before deciding to increase the simulation length.\n", "\n", "## When not to increase simulation length\n", "\n", @@ -154,7 +154,7 @@ "\n", "## Changing the production length\n", "\n", - "The following example shows how to change the production length for a `HybridTopologyProtocol` object. In this example, we will double the production length from `5 ns` to `10 ns` per lambda window. Note the use of units from the `openff.units` package to specify the value is in nanoseconds. We will then calculate the new total simulation time and compare it to the default settings.\n" + "The following example shows how to change the production length for a `RelativeHybridTopologyProtocol` object. In this example, we will double the production length from `5 ns` to `10 ns` per lambda window. Note the use of units from the `openff.units` package to specify the value is in nanoseconds. We will then calculate the new total simulation time and compare it to the default settings.\n" ], "id": "3ab10c2ecebbd3df" }, @@ -225,15 +225,15 @@ "\n", "## Changing the number of lambda windows\n", "\n", - "The following example shows how to change the number of lambda windows for a `HybridTopologyProtocol` object. In this example, we will double the number of lambda windows from `11` to `22`. We will then calculate the new total simulation time and compare it to the default settings. Note the `lambda_windows` setting only accepts integer values and generates a uniform spacing of lambda values between `0` and `1`.\n" + "The following example shows how to change the number of lambda windows for a `RelativeHybridTopologyProtocol` object. In this example, we will double the number of lambda windows from `11` to `22`. We will then calculate the new total simulation time and compare it to the default settings. Note the `lambda_windows` setting only accepts integer values and generates a uniform spacing of lambda values between `0` and `1`.\n" ], "id": "3d6bb208a91ce6d0" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T10:46:24.818876Z", - "start_time": "2026-10-22T10:46:24.776710Z" + "end_time": "2026-10-22T13:14:05.953017Z", + "start_time": "2026-10-22T13:14:05.920141Z" } }, "cell_type": "code", @@ -242,6 +242,8 @@ "\n", "# Increase the number of lambda windows to 22\n", "more_window_settings.lambda_settings.lambda_windows = 22\n", + "# The number of replicas should be increased to match the number of lambda windows\n", + "more_window_settings.simulation_settings.n_replicas = 22\n", "\n", "# Create a new protocol object with the updated settings\n", "protocol = RelativeHybridTopologyProtocol(settings=more_window_settings)\n", @@ -250,7 +252,7 @@ "\n", "# Calculate the new total simulation time\n", "new_total_sim_time = n_repeats * protocol.settings.lambda_settings.lambda_windows * (protocol.settings.simulation_settings.production_length + protocol.settings.simulation_settings.equilibration_length)\n", - "print(f\"New total simulation time: {new_total_sim_time} ns\")\n", + "print(f\"New total simulation time: {new_total_sim_time.m} ns\")\n", "\n", "# Print the percentage increase in total simulation time\n", "percentage_increase = ((new_total_sim_time - total_sim_time) / total_sim_time) * 100\n", @@ -263,12 +265,20 @@ "output_type": "stream", "text": [ "Updated number of lambda windows: 22\n", - "New total simulation time: 396.0 nanosecond ns\n", + "New total simulation time: 396.0 ns\n", "Percentage increase in total simulation time: 100.00%\n" ] } ], - "execution_count": 28 + "execution_count": 37 + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": "", + "id": "2be1cfe17af0153" } ], "metadata": { From ff6d4f2fb26ccc3fa4e0268f358e2a516d8095a5 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Mon, 26 Oct 2026 10:59:36 +0000 Subject: [PATCH 06/11] other protocol details and PR feedback --- cookbook/exploring_rbfe_settings.ipynb | 641 ++++++++++++++++++++++--- 1 file changed, 585 insertions(+), 56 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index 767d5b4..3942ff7 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -4,24 +4,33 @@ "metadata": {}, "cell_type": "markdown", "source": [ - "# Exploring HybridTopologyProtocol settings\n", + "# Exploring Free Energy Protocol settings\n", "\n", + "The free energy protocols implemented in **openfe** expose a number of simulation settings that can affect the convergence properties, accuracy and precision of an alchemical transformation. In this notebook, we will explore changes to common settings, advising when you might want to consider changing them and how to change them in practice. For a detailed explanation of the individual protocols and their limitations see the [protocol documentation](https://docs.openfree.energy/en/latest/guide/protocols/) and references therein.\n" + ], + "id": "a3f965b58be21979" + }, + { + "metadata": {}, + "cell_type": "markdown", + "id": "a1b2c3d4e5f60010", + "source": [ + "## Sampling and performance settings\n", "\n", - "The `RelativeHybridTopologyProtocol` exposes a number of simulation settings that can affect the convergence properties, accuracy and precision of an alchemical transformation. In this notebook, we will explore changes to common settings, advising when you might want to consider changing them and how to change them in practice. For a detailed explanation of the protocol and its limitations see the [protocol documentation](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#relative-hybrid-topology-protocol) and references therein.\n", + "This section covers the settings that control how much sampling is performed and therefore the computational cost and precision of a calculation.\n", "\n", - "## Inspecting the default settings\n", + "### Inspecting the default settings\n", "\n", - "The default settings for the `RelativeHybridTopologyProtocol` are designed to be a good starting point for most alchemical transformations with the goal of obtaining reliable estimates with sufficient precision for decision-making at a reasonable computational cost. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", + "The default settings for each of the protocols are designed to be a good starting point for most alchemical transformations with the goal of obtaining reliable estimates with sufficient precision for decision-making at a reasonable computational cost. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", "\n", - "The default settings can be inspected by:" - ], - "id": "a3f965b58be21979" + "The default settings for the `RelativeHybridTopologyProtocol` can be inspected by:" + ] }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T10:46:22.362791Z", - "start_time": "2026-10-22T10:46:22.275865Z" + "end_time": "2026-10-23T14:46:33.622399Z", + "start_time": "2026-10-23T14:46:33.546091Z" } }, "cell_type": "code", @@ -32,20 +41,26 @@ ], "id": "5690f07a5cd27245", "outputs": [], - "execution_count": 24 + "execution_count": 15 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "For a complete explanation of each setting and how to customize the object see [Choose and Configure a Protocol](https://docs.openfree.energy/en/latest/cookbook/choose_protocol.html).\n", + "
\n", + "For a complete explanation of each setting and how to customize the object see Choose and Configure a Protocol.\n", + "
\n", + "\n", "\n", - "The main parameters that determine the amount of sampling performed and the associated cost for the hybrid topology protocol are:\n", + "The main parameters that determine the amount of sampling performed and the associated cost for these protocols are:\n", "\n", - "- **Lambda windows**: the number of discrete alchemical states simulated during the transformation.\n", - "- **Production length**: the length of the production simulation performed at each lambda window.\n", - "- **Equilibration length**: the length of the equilibration simulation performed at each lambda window.\n", - "- **Protocol repeats**: the number of independent executions of the protocol.\n", + "| Setting | Description | `RelativeHybridTopologyProtocol` | `AbsoluteBindingProtocol` and `SepTopProtocol` |\n", + "|---------------------------------|--------------------------------------------------------------------------------------------------------------------------|--------------------------------------------|--------------------------------------------------------------|\n", + "| **Lambda windows** | The number of discrete alchemical states simulated during the transformation. | `lambda_settings.lambda_windows` | `{complex,solvent}_lambda_settings` |\n", + "| **Production length** | The length of the production simulation performed at each lambda window. | `simulation_settings.production_length` | `{complex,solvent}_simulation_settings.production_length` |\n", + "| **Equilibration length** | The length of the equilibration simulation performed at each lambda window using the selected multisate sampling method. | `simulation_settings.equilibration_length` | `{complex,solvent}_simulation_settings.equilibration_length` |\n", + "| **Pre-equilibration MD length** | The length of the single MD simulation per leg used to determine the restraint atoms in complex legs. | Not applicable | `{complex,solvent}_equil_simulation_settings` |\n", + "| **Protocol repeats** | The number of independent executions of the protocol. | `protocol_repeats` | `protocol_repeats` |\n", "\n", "The appropriate changes to these settings are transformation dependent with no single setting guaranteeing sufficient convergence or accuracy in general. These settings provide a good starting point to investigate the affect of additional sampling on the transformation under study however, but will also increase the computational cost of the simulation. The following example retrieves them from the default settings:" ], @@ -54,8 +69,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T13:12:34.530752Z", - "start_time": "2026-10-22T13:12:34.485473Z" + "end_time": "2026-10-23T14:49:05.032830Z", + "start_time": "2026-10-23T14:49:04.965809Z" } }, "cell_type": "code", @@ -83,13 +98,125 @@ ] } ], - "execution_count": 35 + "execution_count": 16 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "## Estimating the computational cost of a protocol\n", + "#### Other protocols\n", + "\n", + "Similar settings exist for the `AbsoluteBindingProtocol` and `SepTopProtocol`, but they are stored per leg (`complex_*` and `solvent_*`) as those protocols execute both `complex` and `solvent` legs within a single transformation. In addition, each leg has a `*_equil_simulation_settings` object controlling an initial MD simulation used to determine the restraint atoms in the complex legs, which is described in the cost estimate below.\n", + "\n", + "**AbsoluteBindingProtocol**" + ], + "id": "a1b2c3d4e5f60001" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T09:57:16.283825Z", + "start_time": "2026-10-26T09:57:15.875517Z" + } + }, + "cell_type": "code", + "source": [ + "from openfe.protocols.openmm_afe import AbsoluteBindingProtocol\n", + "\n", + "abfe_default_settings = AbsoluteBindingProtocol.default_settings()\n", + "\n", + "print(f\"Independent repeats: {abfe_default_settings.protocol_repeats}\")\n", + "for leg in [\"complex\", \"solvent\"]:\n", + " equil = getattr(abfe_default_settings, f\"{leg}_equil_simulation_settings\")\n", + " sampling = getattr(abfe_default_settings, f\"{leg}_simulation_settings\")\n", + "\n", + " print(f\"{leg}:\")\n", + " print(f\" Lambda windows: {sampling.n_replicas}\")\n", + " print(f\" Production length per state: {sampling.production_length.m} ns\")\n", + " print(f\" Equilibration length per state: {sampling.equilibration_length.m} ns\")\n", + " print(f\" Pre-equilibration MD length: {(equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length).m} ns\")" + ], + "id": "a1b2c3d4e5f60002", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Independent repeats: 3\n", + "complex:\n", + " Lambda windows: 30\n", + " Production length per state: 10.0 ns\n", + " Equilibration length per state: 1 ns\n", + " Pre-equilibration MD length: 5.75 ns\n", + "solvent:\n", + " Lambda windows: 14\n", + " Production length per state: 10.0 ns\n", + " Equilibration length per state: 1.0 ns\n", + " Pre-equilibration MD length: 0.8 ns\n" + ] + } + ], + "execution_count": 49 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "**SepTopProtocol**" + ], + "id": "a1b2c3d4e5f60003" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T09:57:39.281431Z", + "start_time": "2026-10-26T09:57:39.219301Z" + } + }, + "cell_type": "code", + "source": [ + "from openfe.protocols.openmm_septop import SepTopProtocol\n", + "\n", + "septop_default_settings = SepTopProtocol.default_settings()\n", + "\n", + "print(f\"Independent repeats: {septop_default_settings.protocol_repeats}\")\n", + "for leg in [\"complex\", \"solvent\"]:\n", + " equil = getattr(septop_default_settings, f\"{leg}_equil_simulation_settings\")\n", + " sampling = getattr(septop_default_settings, f\"{leg}_simulation_settings\")\n", + "\n", + " print(f\"{leg}:\")\n", + " print(f\" Lambda windows: {sampling.n_replicas}\")\n", + " print(f\" Production length per state: {sampling.production_length.m} ns\")\n", + " print(f\" Equilibration length per state: {sampling.equilibration_length.m} ns\")\n", + " print(f\" Pre-equilibration MD length: {(equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length).m} ns\")" + ], + "id": "a1b2c3d4e5f60004", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Independent repeats: 3\n", + "complex:\n", + " Lambda windows: 19\n", + " Production length per state: 10.0 ns\n", + " Equilibration length per state: 1.0 ns\n", + " Pre-equilibration MD length: 2.2 ns\n", + "solvent:\n", + " Lambda windows: 27\n", + " Production length per state: 10.0 ns\n", + " Equilibration length per state: 1.0 ns\n", + " Pre-equilibration MD length: 2.2 ns\n" + ] + } + ], + "execution_count": 50 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Estimating the computational cost of a protocol\n", "\n", "A useful first approximation of the computational cost of a protocol execution is the total simulation time in nanoseconds across all lambda windows and repeats. This can be calculated as:\n", "\n", @@ -102,14 +229,22 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T11:20:41.651615Z", - "start_time": "2026-10-22T11:20:41.533380Z" + "end_time": "2026-10-26T09:57:52.628291Z", + "start_time": "2026-10-26T09:57:52.573302Z" } }, "cell_type": "code", "source": [ - "total_sim_time = n_repeats * n_windows * (production_length + equilibration_length)\n", - "print(f\"Total simulation time: {total_sim_time.m} ns\")" + "def estimate_rbfe_cost(settings):\n", + " n_repeats = settings.protocol_repeats\n", + " n_windows = settings.lambda_settings.lambda_windows\n", + " production_length = settings.simulation_settings.production_length\n", + " equilibration_length = settings.simulation_settings.equilibration_length\n", + "\n", + " total_sim_time = n_repeats * n_windows * (production_length + equilibration_length)\n", + " return total_sim_time\n", + "total_sim_time_per_leg = estimate_rbfe_cost(settings)\n", + "print(f\"Total simulation time: {total_sim_time_per_leg.m} ns\")" ], "id": "fd81e13e2c34684d", "outputs": [ @@ -121,48 +256,172 @@ ] } ], - "execution_count": 29 + "execution_count": 51 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "The `RelativeHybridTopologyProtocol` simulates the complex and solvent legs of a transformation independently, so if both use the same settings the total cost of the protocol is double the above estimate. It is useful to consider this estimate when deciding how to adjust the settings for a given transformation.\n", + "The `RelativeHybridTopologyProtocol` simulates the complex and solvent legs of a transformation independently, so if both use the same settings the total cost of the protocol is double the above estimate. It is useful to consider this estimate when deciding how to adjust the settings for a given transformation." + ], + "id": "3ab10c2ecebbd3df" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Estimating the cost of the ABFE and SepTop protocols\n", "\n", - "# Increasing simulation length\n", + "The `AbsoluteBindingProtocol` and `SepTopProtocol` have a slightly different cost profile. Before the production simulation, each leg first runs a single, longer MD equilibration (NVT equilibration, NPT equilibration and a production run) that is used to determine the restraint atoms for the complex leg and seed the solvent leg. This simulation is run once per leg and repeat, **not** once per lambda window, so the total simulation time becomes:\n", + "\n", + "$$T_{total} = N_{repeats} \\times \\sum_{legs} \\left[ T_{MD} + N_{\\lambda windows} \\times (T_{production} + T_{equilibration}) \\right]$$\n", + "\n", + "where $T_{MD}$ is the sum of the NVT equilibration, equilibration and production lengths of the leg's `*_equil_simulation_settings`, and the remaining terms come from its `*_simulation_settings`. Both protocols require the number of replicas to match the number of lambda windows, so `n_replicas` gives the number of windows. The legs are `complex` and `solvent`." + ], + "id": "5017f2cf" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "**AbsoluteBindingProtocol**" + ], + "id": "50aaaf76" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T10:26:00.830205Z", + "start_time": "2026-10-26T10:26:00.795098Z" + } + }, + "cell_type": "code", + "source": [ + "from openff.units import unit\n", + "\n", + "def estimate_abfe_cost(settings):\n", + " abfe_total_sim_time = 0 * unit.nanoseconds\n", + " for leg in [\"complex\", \"solvent\"]:\n", + " equil = getattr(settings, f\"{leg}_equil_simulation_settings\")\n", + " sampling = getattr(settings, f\"{leg}_simulation_settings\")\n", + "\n", + " md_time = equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length\n", + " replica_exchange_time = sampling.n_replicas * (sampling.production_length + sampling.equilibration_length)\n", + " leg_time = md_time + replica_exchange_time\n", + " print(f\"{leg}: {md_time.m} ns MD + {replica_exchange_time.m} ns replica exchange = {leg_time.m} ns\")\n", + " abfe_total_sim_time += leg_time\n", + "\n", + " abfe_total_sim_time *= settings.protocol_repeats\n", + " print(f\"Total simulation time: {abfe_total_sim_time.m:.2f} ns\")\n", + " return abfe_total_sim_time\n", + "\n", + "default_abfe_cost = estimate_abfe_cost(abfe_default_settings)" + ], + "id": "0f8fb861", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "complex: 5.75 ns MD + 330.0 ns replica exchange = 335.75 ns\n", + "solvent: 0.8 ns MD + 154.0 ns replica exchange = 154.8 ns\n", + "Total simulation time: 1471.65 ns\n" + ] + } + ], + "execution_count": 81 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "**SepTopProtocol**" + ], + "id": "b004911d" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T10:23:53.381934Z", + "start_time": "2026-10-26T10:23:53.342287Z" + } + }, + "cell_type": "code", + "source": [ + "def estimate_septop_cost(settings):\n", + " septop_total_sim_time = 0 * unit.nanoseconds\n", + " for leg in [\"complex\", \"solvent\"]:\n", + " equil = getattr(settings, f\"{leg}_equil_simulation_settings\")\n", + " sampling = getattr(settings, f\"{leg}_simulation_settings\")\n", + "\n", + " md_time = equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length\n", + " replica_exchange_time = sampling.n_replicas * (sampling.production_length + sampling.equilibration_length)\n", + " leg_time = md_time + replica_exchange_time\n", + " print(f\"{leg}: {md_time.m} ns MD + {replica_exchange_time.m} ns replica exchange = {leg_time.m} ns\")\n", + " septop_total_sim_time += leg_time\n", + "\n", + " septop_total_sim_time *= settings.protocol_repeats\n", + " print(f\"Total simulation time: {septop_total_sim_time.m:.2f} ns\")\n", + " return septop_total_sim_time\n", + "\n", + "default_septop_cost = estimate_septop_cost(septop_default_settings)" + ], + "id": "85ed4076", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "complex: 2.2 ns MD + 209.0 ns replica exchange = 211.2 ns\n", + "solvent: 2.2 ns MD + 297.0 ns replica exchange = 299.2 ns\n", + "Total simulation time: 1531.20 ns\n" + ] + } + ], + "execution_count": 74 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Increasing simulation length\n", "\n", "This section describes how to increase the amount of production simulation performed at each lambda window and estimate the increase in computational cost.\n", "\n", - "## When to increase simulation length\n", + "
\n", + "Warning: Extending an existing simulation is not currently supported. Increasing the simulation length after running with default settings means the protocol must be run again from the start with the new settings, so the full cost of the longer simulation is incurred.\n", + "
\n", + "\n", + "#### When to increase simulation length\n", "\n", "Consider increasing the simulation length **before** running the simulation when:\n", "\n", - "- The transformation involves large structural changes due to for example changes in formal charge, you may also need to consider increasing the number of lambda windows see the next section for more information.\n", + "- The transformation involves large structural changes due to for example changes in formal charge. In this case you may also need to consider increasing the number of lambda windows see the next section for more information.\n", "\n", "Consider increasing the simulation length **after** running the simulation when:\n", "\n", "- Free-energy estimates differ substantially between independent repeats of the protocol.\n", - "- Forward and reverse free-energy estimates disagree by more than 1kbT at over 50% of the simulation length and their uncertainties do not overlap.\n", + "- Forward and reverse free-energy estimates disagree by more than 1kbT at over 50% of the simulation length and their uncertainties do not overlap.[[1]](https://livecomsjournal.org/index.php/livecoms/article/view/v2i1e18378)\n", "- A large number of samples are discarded during analysis due to the automated equilibration and decorrelation methods.\n", "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows, you may also need to consider increasing the number of lambda windows see the next section for more information.\n", "\n", - "This information is available through the automated analysis performed by the `RelativeHybridTopologyProtocol` see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for more information, it is recommended to review the analysis metrics before deciding to increase the simulation length.\n", + "This information is available through the automated analysis performed by the protocols the see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for the `RelativeHybridTopologyProtocol` for more information. It is recommended to review the analysis metrics before deciding to increase the simulation length.\n", "\n", - "## When not to increase simulation length\n", + "#### When not to increase simulation length\n", "\n", - "Increasing the simulation length is not always the best solution. For example poor overlap between neighboring lambda windows may require changes to the number of lambda windows or the spacing between them.\n", + "Increasing the simulation length is not always the best solution. For example changes to the number of lambda windows or the spacing between them can sometimes be more effective when you have poor overlap between neighboring lambda windows and in some cases a combination of both maybe required.\n", "\n", - "## Changing the production length\n", + "#### Changing the production length\n", "\n", "The following example shows how to change the production length for a `RelativeHybridTopologyProtocol` object. In this example, we will double the production length from `5 ns` to `10 ns` per lambda window. Note the use of units from the `openff.units` package to specify the value is in nanoseconds. We will then calculate the new total simulation time and compare it to the default settings.\n" ], - "id": "3ab10c2ecebbd3df" + "id": "771ba6b9" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T11:21:00.199353Z", - "start_time": "2026-10-22T11:21:00.150806Z" + "end_time": "2026-10-26T10:23:54.390009Z", + "start_time": "2026-10-26T10:23:54.335215Z" } }, "cell_type": "code", @@ -180,11 +439,11 @@ "print(f\"Updated production length: {protocol.settings.simulation_settings.production_length.m} ns\")\n", "\n", "# Calculate the new total simulation time\n", - "new_total_sim_time = n_repeats * n_windows * (protocol.settings.simulation_settings.production_length + protocol.settings.simulation_settings.equilibration_length)\n", - "print(f\"New total simulation time: {new_total_sim_time.m} ns\")\n", + "new_total_sim_time_per_leg = estimate_rbfe_cost(protocol.settings)\n", + "print(f\"New total simulation time: {new_total_sim_time_per_leg.m} ns\")\n", "\n", "# Print the percentage increase in total simulation time\n", - "percentage_increase = ((new_total_sim_time - total_sim_time) / total_sim_time) * 100\n", + "percentage_increase = ((new_total_sim_time_per_leg - total_sim_time_per_leg) / total_sim_time_per_leg) * 100\n", "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" ], "id": "2cf9b3c7030df36f", @@ -199,17 +458,130 @@ ] } ], - "execution_count": 31 + "execution_count": 75 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "# Increasing the number of lambda windows\n", + "##### Other protocols\n", + "\n", + "The same idea applies to each of the other protocols, only the location of the setting changes. Protocols that simulate several legs (e.g. complex and solvent) store the settings per leg, so you must change each leg to be extended individually.\n", + "\n", + "**AbsoluteBindingProtocol**" + ], + "id": "d4f1a7c2e9b34a60" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T10:23:55.747906Z", + "start_time": "2026-10-26T10:23:55.703515Z" + } + }, + "cell_type": "code", + "source": [ + "from openfe.protocols.openmm_afe import AbsoluteBindingProtocol\n", + "\n", + "abfe_long_settings = AbsoluteBindingProtocol.default_settings()\n", + "\n", + "# Settings are per leg, so extend both legs\n", + "abfe_long_settings.complex_simulation_settings.production_length = 20 * unit.nanoseconds\n", + "abfe_long_settings.solvent_simulation_settings.production_length = 20 * unit.nanoseconds\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "abfe_protocol = AbsoluteBindingProtocol(settings=abfe_long_settings)\n", + "\n", + "print(f\"Updated complex production length: {abfe_protocol.settings.complex_simulation_settings.production_length.m} ns\")\n", + "print(f\"Updated solvent production length: {abfe_protocol.settings.solvent_simulation_settings.production_length.m} ns\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_abfe_cost = estimate_abfe_cost(abfe_long_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_abfe_cost - default_abfe_cost) / default_abfe_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "e1a2b3c4d5f60718", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated complex production length: 20 ns\n", + "Updated solvent production length: 20 ns\n", + "complex: 5.75 ns MD + 630 ns replica exchange = 635.75 ns\n", + "solvent: 0.8 ns MD + 294.0 ns replica exchange = 294.8 ns\n", + "Total simulation time: 2791.65 ns\n", + "Percentage increase in total simulation time: 89.70%\n" + ] + } + ], + "execution_count": 76 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "**SepTopProtocol**" + ], + "id": "f2b3c4d5e6a71829" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T10:23:57.529452Z", + "start_time": "2026-10-26T10:23:57.482671Z" + } + }, + "cell_type": "code", + "source": [ + "from openfe.protocols.openmm_septop import SepTopProtocol\n", + "\n", + "septop_long_settings = SepTopProtocol.default_settings()\n", + "\n", + "septop_long_settings.complex_simulation_settings.production_length = 20 * unit.nanoseconds\n", + "septop_long_settings.solvent_simulation_settings.production_length = 20 * unit.nanoseconds\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "septop_protocol = SepTopProtocol(settings=septop_long_settings)\n", + "\n", + "print(f\"Updated complex production length: {septop_protocol.settings.complex_simulation_settings.production_length.m} ns\")\n", + "print(f\"Updated solvent production length: {septop_protocol.settings.solvent_simulation_settings.production_length.m} ns\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_septop_cost = estimate_septop_cost(septop_long_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_septop_cost - default_septop_cost) / default_septop_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "a3c4d5e6f7b82930", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated complex production length: 20 ns\n", + "Updated solvent production length: 20 ns\n", + "complex: 2.2 ns MD + 399.0 ns replica exchange = 401.2 ns\n", + "solvent: 2.2 ns MD + 567.0 ns replica exchange = 569.2 ns\n", + "Total simulation time: 2911.20 ns\n", + "Percentage increase in total simulation time: 90.13%\n" + ] + } + ], + "execution_count": 77 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Increasing the number of lambda windows\n", "\n", "Increasing the number of lambda windows introduces additional intermediate alchemical states along the transformation which can make the pathway easier to traverse and improve the overlap between neighboring states. This section describes how to increase the number of lambda windows and estimate the increase in computational cost.\n", "\n", - "## When to increase the number of lambda windows\n", + "#### When to increase the number of lambda windows\n", "\n", "Consider increasing the number of lambda windows **before** running the simulation when:\n", "\n", @@ -219,11 +591,11 @@ "\n", "- Neighboring lambda windows have poor phase-space overlap or replica mixing statistics.\n", "\n", - "## When not to increase the number of lambda windows\n", + "#### When not to increase the number of lambda windows\n", "\n", - "Increasing the number of lambda windows is not always the best solution, slow changes in ligand conformation or binding site rearrangements may require longer simulation lengths to achieve convergence. In these cases, increasing the number of lambda windows may not improve convergence and will increase the computational cost of the protocol.\n", + "Increasing the number of lambda windows is not always the best solution, slow changes in ligand conformation or binding site rearrangements may require longer simulation lengths to achieve convergence.\n", "\n", - "## Changing the number of lambda windows\n", + "#### Changing the number of lambda windows\n", "\n", "The following example shows how to change the number of lambda windows for a `RelativeHybridTopologyProtocol` object. In this example, we will double the number of lambda windows from `11` to `22`. We will then calculate the new total simulation time and compare it to the default settings. Note the `lambda_windows` setting only accepts integer values and generates a uniform spacing of lambda values between `0` and `1`.\n" ], @@ -232,8 +604,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-22T13:14:05.953017Z", - "start_time": "2026-10-22T13:14:05.920141Z" + "end_time": "2026-10-26T10:24:06.008414Z", + "start_time": "2026-10-26T10:24:05.974089Z" } }, "cell_type": "code", @@ -242,7 +614,7 @@ "\n", "# Increase the number of lambda windows to 22\n", "more_window_settings.lambda_settings.lambda_windows = 22\n", - "# The number of replicas should be increased to match the number of lambda windows\n", + "# The number of replicas must be increased to match the number of lambda windows\n", "more_window_settings.simulation_settings.n_replicas = 22\n", "\n", "# Create a new protocol object with the updated settings\n", @@ -251,11 +623,11 @@ "print(f\"Updated number of lambda windows: {protocol.settings.lambda_settings.lambda_windows}\")\n", "\n", "# Calculate the new total simulation time\n", - "new_total_sim_time = n_repeats * protocol.settings.lambda_settings.lambda_windows * (protocol.settings.simulation_settings.production_length + protocol.settings.simulation_settings.equilibration_length)\n", - "print(f\"New total simulation time: {new_total_sim_time.m} ns\")\n", + "new_total_sim_time_per_leg = estimate_rbfe_cost(protocol.settings)\n", + "print(f\"New total simulation time: {new_total_sim_time_per_leg.m} ns\")\n", "\n", "# Print the percentage increase in total simulation time\n", - "percentage_increase = ((new_total_sim_time - total_sim_time) / total_sim_time) * 100\n", + "percentage_increase = ((new_total_sim_time_per_leg - total_sim_time_per_leg) / total_sim_time_per_leg) * 100\n", "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" ], "id": "a8a0a2b5ba644eb4", @@ -270,15 +642,172 @@ ] } ], - "execution_count": 37 + "execution_count": 78 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "##### Other protocols\n", + "\n", + "The `AbsoluteBindingProtocol` and `SepTopProtocol` allow for exact control of the lambda schedule of each leg requiring an explicitly defined list of values, one entry per window. All of the lists for a leg must have the same length, and `n_replicas` in the matching `*_simulation_settings` must be set to the number of windows. The examples below change the `complex` leg, the `solvent` leg can be changed in the same way using `solvent_lambda_settings` and `solvent_simulation_settings`.\n", + "\n", + "**AbsoluteBindingProtocol**\n", + "\n", + "The schedule is defined by `lambda_elec`, `lambda_vdw` and `lambda_restraints`." + ], + "id": "b1c2d3e4f5a60001" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T10:27:37.430923Z", + "start_time": "2026-10-26T10:27:37.389418Z" + } + }, + "cell_type": "code", + "source": [ + "abfe_window_settings = AbsoluteBindingProtocol.default_settings()\n", + "\n", + "# Define the new lambda schedule for the complex leg by adding 3 more charge windows\n", + "abfe_window_settings.complex_lambda_settings.lambda_elec = [0.0] * 6 + [0.1, 0.2, 0.3, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.85, 0.9] + [1.0] * 15\n", + "# Extra vdW windows are added to match the length of the new charge schedule\n", + "abfe_window_settings.complex_lambda_settings.lambda_vdw = [0.0] * 19 + [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]\n", + "# Extra restraint windows are added to match the length of the new charge schedule\n", + "abfe_window_settings.complex_lambda_settings.lambda_restraints = [0.0, 0.2, 0.4, 0.6, 0.8] + [1.0] * 28\n", + "\n", + "# The number of replicas must match the number of lambda windows\n", + "abfe_window_settings.complex_simulation_settings.n_replicas = len(abfe_window_settings.complex_lambda_settings.lambda_vdw)\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "abfe_protocol = AbsoluteBindingProtocol(settings=abfe_window_settings)\n", + "\n", + "# Calculate the new total simulation time\n", + "new_abfe_cost = estimate_abfe_cost(abfe_window_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_abfe_cost - default_abfe_cost) / default_abfe_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "b1c2d3e4f5a60002", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "complex: 5.75 ns MD + 363.0 ns replica exchange = 368.75 ns\n", + "solvent: 0.8 ns MD + 154.0 ns replica exchange = 154.8 ns\n", + "Total simulation time: 1570.65 ns\n", + "Percentage increase in total simulation time: 6.73%\n" + ] + } + ], + "execution_count": 82 }, { "metadata": {}, + "cell_type": "markdown", + "source": [ + "
\n", + "Tip: Changing the number of lambda windows or the lambda schedule for the AbsoluteBindingProtocol and SepTopProtocol by hand can be error prone. The protocols provide a validation method to check that the lambda schedules are valid for the leg in question which can be used to check your changes before running the protocol. See the example below.\n", + "
" + ], + "id": "2290219ca360bbd3" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T10:57:39.961356Z", + "start_time": "2026-10-26T10:57:39.949172Z" + } + }, + "cell_type": "code", + "source": [ + "# Validate the lambda schedules for the complex leg\n", + "AbsoluteBindingProtocol._validate_lambda_schedule(\n", + " abfe_window_settings.complex_lambda_settings,\n", + " abfe_window_settings.complex_simulation_settings\n", + ")" + ], + "id": "6c9de9ae4d246d9c", + "outputs": [], + "execution_count": 88 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "**SepTopProtocol**\n", + "\n", + "The schedule is defined by the `A` and `B` states of each component: `lambda_elec_A`, `lambda_elec_B`, `lambda_vdw_A`, `lambda_vdw_B`, `lambda_restraints_A` and `lambda_restraints_B`. Like the `AbsoluteBindingProtocol` all the lambda value lists for a leg must be the same length." + ], + "id": "b1c2d3e4f5a60003" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T10:58:15.164424Z", + "start_time": "2026-10-26T10:58:15.108837Z" + } + }, + "cell_type": "code", + "source": [ + "septop_window_settings = SepTopProtocol.default_settings()\n", + "\n", + "# Define the new lambda schedule for the complex leg adding 3 more charge windows for A and B\n", + "septop_window_settings.complex_lambda_settings.lambda_elec_A = [0.0] * 8 + [0.25, 0.4, 0.5, 0.6, 0.75, 0.8] + [1.0] * 8\n", + "septop_window_settings.complex_lambda_settings.lambda_elec_B = [1.0] * 8 + [0.8, 0.75, 0.6, 0.5, 0.4, 0.25] + [0.0] * 8\n", + "septop_window_settings.complex_lambda_settings.lambda_vdw_A = [0.0] * 15 + [0.14, 0.28, 0.42, 0.57, 0.71, 0.85, 1.0]\n", + "septop_window_settings.complex_lambda_settings.lambda_vdw_B = [1.0, 0.85, 0.71, 0.57, 0.42, 0.28, 0.14] + [0.0] * 15\n", + "septop_window_settings.complex_lambda_settings.lambda_restraints_A = [0.0, 0.05, 0.1, 0.3, 0.5, 0.75] + [1.0] * 16\n", + "septop_window_settings.complex_lambda_settings.lambda_restraints_B = [1.0] * 16 + [0.75, 0.5, 0.3, 0.1, 0.05, 0.0]\n", + "\n", + "# The number of replicas must match the number of lambda windows\n", + "septop_window_settings.complex_simulation_settings.n_replicas = len(septop_window_settings.complex_lambda_settings.lambda_vdw_A)\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "septop_protocol = SepTopProtocol(settings=septop_window_settings)\n", + "\n", + "# Calculate the new total simulation time\n", + "new_septop_cost = estimate_septop_cost(septop_window_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_septop_cost - default_septop_cost) / default_septop_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "b1c2d3e4f5a60004", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "complex: 2.2 ns MD + 242.0 ns replica exchange = 244.2 ns\n", + "solvent: 2.2 ns MD + 297.0 ns replica exchange = 299.2 ns\n", + "Total simulation time: 1630.20 ns\n", + "Percentage increase in total simulation time: 6.47%\n" + ] + } + ], + "execution_count": 90 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T10:59:08.894162Z", + "start_time": "2026-10-26T10:59:08.880003Z" + } + }, "cell_type": "code", + "source": [ + "# Validate the lambda schedules for the complex leg\n", + "SepTopProtocol._validate_lambda_schedule(\n", + " septop_window_settings.complex_lambda_settings,\n", + " septop_window_settings.complex_simulation_settings\n", + ")" + ], + "id": "32e3f17b3a1e849", "outputs": [], - "execution_count": null, - "source": "", - "id": "2be1cfe17af0153" + "execution_count": 91 } ], "metadata": { From f4207db6fd50f6597aa6fa3aaa1564af7d0ad9c8 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Mon, 26 Oct 2026 11:45:26 +0000 Subject: [PATCH 07/11] more cleanup --- cookbook/exploring_rbfe_settings.ipynb | 141 ++++++++++++------------- 1 file changed, 69 insertions(+), 72 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index 3942ff7..80c307c 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -29,8 +29,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-23T14:46:33.622399Z", - "start_time": "2026-10-23T14:46:33.546091Z" + "end_time": "2026-10-26T11:45:14.335951Z", + "start_time": "2026-10-26T11:45:01.020777Z" } }, "cell_type": "code", @@ -41,7 +41,7 @@ ], "id": "5690f07a5cd27245", "outputs": [], - "execution_count": 15 + "execution_count": 1 }, { "metadata": {}, @@ -69,8 +69,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-23T14:49:05.032830Z", - "start_time": "2026-10-23T14:49:04.965809Z" + "end_time": "2026-10-26T11:45:14.395446Z", + "start_time": "2026-10-26T11:45:14.352357Z" } }, "cell_type": "code", @@ -98,7 +98,7 @@ ] } ], - "execution_count": 16 + "execution_count": 2 }, { "metadata": {}, @@ -115,8 +115,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T09:57:16.283825Z", - "start_time": "2026-10-26T09:57:15.875517Z" + "end_time": "2026-10-26T11:45:16.502740Z", + "start_time": "2026-10-26T11:45:14.398886Z" } }, "cell_type": "code", @@ -154,9 +154,24 @@ " Equilibration length per state: 1.0 ns\n", " Pre-equilibration MD length: 0.8 ns\n" ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/joshua/mambaforge/envs/openfe_dev/lib/python3.12/site-packages/Bio/Application/__init__.py:39: BiopythonDeprecationWarning: The Bio.Application modules and modules relying on it have been deprecated.\n", + "\n", + "Due to the on going maintenance burden of keeping command line application\n", + "wrappers up to date, we have decided to deprecate and eventually remove these\n", + "modules.\n", + "\n", + "We instead now recommend building your command line and invoking it directly\n", + "with the subprocess module.\n", + " warnings.warn(\n" + ] } ], - "execution_count": 49 + "execution_count": 3 }, { "metadata": {}, @@ -169,8 +184,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T09:57:39.281431Z", - "start_time": "2026-10-26T09:57:39.219301Z" + "end_time": "2026-10-26T11:45:16.632171Z", + "start_time": "2026-10-26T11:45:16.525149Z" } }, "cell_type": "code", @@ -210,7 +225,7 @@ ] } ], - "execution_count": 50 + "execution_count": 4 }, { "metadata": {}, @@ -229,8 +244,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T09:57:52.628291Z", - "start_time": "2026-10-26T09:57:52.573302Z" + "end_time": "2026-10-26T11:45:16.681723Z", + "start_time": "2026-10-26T11:45:16.656346Z" } }, "cell_type": "code", @@ -256,7 +271,7 @@ ] } ], - "execution_count": 51 + "execution_count": 5 }, { "metadata": {}, @@ -291,15 +306,15 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:26:00.830205Z", - "start_time": "2026-10-26T10:26:00.795098Z" + "end_time": "2026-10-26T11:45:16.787042Z", + "start_time": "2026-10-26T11:45:16.711856Z" } }, "cell_type": "code", "source": [ "from openff.units import unit\n", "\n", - "def estimate_abfe_cost(settings):\n", + "def estimate_abfe_septop_cost(settings):\n", " abfe_total_sim_time = 0 * unit.nanoseconds\n", " for leg in [\"complex\", \"solvent\"]:\n", " equil = getattr(settings, f\"{leg}_equil_simulation_settings\")\n", @@ -315,7 +330,7 @@ " print(f\"Total simulation time: {abfe_total_sim_time.m:.2f} ns\")\n", " return abfe_total_sim_time\n", "\n", - "default_abfe_cost = estimate_abfe_cost(abfe_default_settings)" + "default_abfe_cost = estimate_abfe_septop_cost(abfe_default_settings)" ], "id": "0f8fb861", "outputs": [ @@ -329,7 +344,7 @@ ] } ], - "execution_count": 81 + "execution_count": 6 }, { "metadata": {}, @@ -342,30 +357,12 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:23:53.381934Z", - "start_time": "2026-10-26T10:23:53.342287Z" + "end_time": "2026-10-26T11:45:16.850622Z", + "start_time": "2026-10-26T11:45:16.790688Z" } }, "cell_type": "code", - "source": [ - "def estimate_septop_cost(settings):\n", - " septop_total_sim_time = 0 * unit.nanoseconds\n", - " for leg in [\"complex\", \"solvent\"]:\n", - " equil = getattr(settings, f\"{leg}_equil_simulation_settings\")\n", - " sampling = getattr(settings, f\"{leg}_simulation_settings\")\n", - "\n", - " md_time = equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length\n", - " replica_exchange_time = sampling.n_replicas * (sampling.production_length + sampling.equilibration_length)\n", - " leg_time = md_time + replica_exchange_time\n", - " print(f\"{leg}: {md_time.m} ns MD + {replica_exchange_time.m} ns replica exchange = {leg_time.m} ns\")\n", - " septop_total_sim_time += leg_time\n", - "\n", - " septop_total_sim_time *= settings.protocol_repeats\n", - " print(f\"Total simulation time: {septop_total_sim_time.m:.2f} ns\")\n", - " return septop_total_sim_time\n", - "\n", - "default_septop_cost = estimate_septop_cost(septop_default_settings)" - ], + "source": "default_septop_cost = estimate_abfe_septop_cost(septop_default_settings)", "id": "85ed4076", "outputs": [ { @@ -378,7 +375,7 @@ ] } ], - "execution_count": 74 + "execution_count": 7 }, { "metadata": {}, @@ -405,7 +402,7 @@ "- A large number of samples are discarded during analysis due to the automated equilibration and decorrelation methods.\n", "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows, you may also need to consider increasing the number of lambda windows see the next section for more information.\n", "\n", - "This information is available through the automated analysis performed by the protocols the see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for the `RelativeHybridTopologyProtocol` for more information. It is recommended to review the analysis metrics before deciding to increase the simulation length.\n", + "This information is available through the automated analysis performed by the protocols, see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for the `RelativeHybridTopologyProtocol` for more information. It is recommended to review the analysis metrics before deciding to increase the simulation length.\n", "\n", "#### When not to increase simulation length\n", "\n", @@ -420,8 +417,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:23:54.390009Z", - "start_time": "2026-10-26T10:23:54.335215Z" + "end_time": "2026-10-26T11:45:16.910948Z", + "start_time": "2026-10-26T11:45:16.853625Z" } }, "cell_type": "code", @@ -458,7 +455,7 @@ ] } ], - "execution_count": 75 + "execution_count": 8 }, { "metadata": {}, @@ -475,8 +472,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:23:55.747906Z", - "start_time": "2026-10-26T10:23:55.703515Z" + "end_time": "2026-10-26T11:45:16.965069Z", + "start_time": "2026-10-26T11:45:16.913325Z" } }, "cell_type": "code", @@ -496,7 +493,7 @@ "print(f\"Updated solvent production length: {abfe_protocol.settings.solvent_simulation_settings.production_length.m} ns\")\n", "\n", "# Calculate the new total simulation time\n", - "new_abfe_cost = estimate_abfe_cost(abfe_long_settings)\n", + "new_abfe_cost = estimate_abfe_septop_cost(abfe_long_settings)\n", "\n", "# Print the percentage increase in total simulation time\n", "percentage_increase = ((new_abfe_cost - default_abfe_cost) / default_abfe_cost) * 100\n", @@ -517,7 +514,7 @@ ] } ], - "execution_count": 76 + "execution_count": 9 }, { "metadata": {}, @@ -530,8 +527,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:23:57.529452Z", - "start_time": "2026-10-26T10:23:57.482671Z" + "end_time": "2026-10-26T11:45:17.016519Z", + "start_time": "2026-10-26T11:45:16.983332Z" } }, "cell_type": "code", @@ -550,7 +547,7 @@ "print(f\"Updated solvent production length: {septop_protocol.settings.solvent_simulation_settings.production_length.m} ns\")\n", "\n", "# Calculate the new total simulation time\n", - "new_septop_cost = estimate_septop_cost(septop_long_settings)\n", + "new_septop_cost = estimate_abfe_septop_cost(septop_long_settings)\n", "\n", "# Print the percentage increase in total simulation time\n", "percentage_increase = ((new_septop_cost - default_septop_cost) / default_septop_cost) * 100\n", @@ -571,7 +568,7 @@ ] } ], - "execution_count": 77 + "execution_count": 10 }, { "metadata": {}, @@ -604,8 +601,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:24:06.008414Z", - "start_time": "2026-10-26T10:24:05.974089Z" + "end_time": "2026-10-26T11:45:17.161411Z", + "start_time": "2026-10-26T11:45:17.043534Z" } }, "cell_type": "code", @@ -642,7 +639,7 @@ ] } ], - "execution_count": 78 + "execution_count": 11 }, { "metadata": {}, @@ -661,8 +658,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:27:37.430923Z", - "start_time": "2026-10-26T10:27:37.389418Z" + "end_time": "2026-10-26T11:45:17.218775Z", + "start_time": "2026-10-26T11:45:17.183233Z" } }, "cell_type": "code", @@ -683,7 +680,7 @@ "abfe_protocol = AbsoluteBindingProtocol(settings=abfe_window_settings)\n", "\n", "# Calculate the new total simulation time\n", - "new_abfe_cost = estimate_abfe_cost(abfe_window_settings)\n", + "new_abfe_cost = estimate_abfe_septop_cost(abfe_window_settings)\n", "\n", "# Print the percentage increase in total simulation time\n", "percentage_increase = ((new_abfe_cost - default_abfe_cost) / default_abfe_cost) * 100\n", @@ -702,14 +699,14 @@ ] } ], - "execution_count": 82 + "execution_count": 12 }, { "metadata": {}, "cell_type": "markdown", "source": [ "
\n", - "Tip: Changing the number of lambda windows or the lambda schedule for the AbsoluteBindingProtocol and SepTopProtocol by hand can be error prone. The protocols provide a validation method to check that the lambda schedules are valid for the leg in question which can be used to check your changes before running the protocol. See the example below.\n", + "Tip: Changing the number of lambda windows or the lambda schedule for the AbsoluteBindingProtocol and SepTopProtocol by hand can be error prone. The protocols provide a validation method to check that the lambda schedules are valid for the leg in question, use them to validate your changes and catch issues before execution. See the example below.\n", "
" ], "id": "2290219ca360bbd3" @@ -717,8 +714,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:57:39.961356Z", - "start_time": "2026-10-26T10:57:39.949172Z" + "end_time": "2026-10-26T11:45:17.251783Z", + "start_time": "2026-10-26T11:45:17.234232Z" } }, "cell_type": "code", @@ -731,7 +728,7 @@ ], "id": "6c9de9ae4d246d9c", "outputs": [], - "execution_count": 88 + "execution_count": 13 }, { "metadata": {}, @@ -746,8 +743,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:58:15.164424Z", - "start_time": "2026-10-26T10:58:15.108837Z" + "end_time": "2026-10-26T11:45:17.286028Z", + "start_time": "2026-10-26T11:45:17.264318Z" } }, "cell_type": "code", @@ -769,7 +766,7 @@ "septop_protocol = SepTopProtocol(settings=septop_window_settings)\n", "\n", "# Calculate the new total simulation time\n", - "new_septop_cost = estimate_septop_cost(septop_window_settings)\n", + "new_septop_cost = estimate_abfe_septop_cost(septop_window_settings)\n", "\n", "# Print the percentage increase in total simulation time\n", "percentage_increase = ((new_septop_cost - default_septop_cost) / default_septop_cost) * 100\n", @@ -788,13 +785,13 @@ ] } ], - "execution_count": 90 + "execution_count": 14 }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T10:59:08.894162Z", - "start_time": "2026-10-26T10:59:08.880003Z" + "end_time": "2026-10-26T11:45:17.315038Z", + "start_time": "2026-10-26T11:45:17.307686Z" } }, "cell_type": "code", @@ -807,7 +804,7 @@ ], "id": "32e3f17b3a1e849", "outputs": [], - "execution_count": 91 + "execution_count": 15 } ], "metadata": { From 70308bca13d47a1dae81dfca37a1f8e90a4dceb1 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Mon, 26 Oct 2026 13:11:32 +0000 Subject: [PATCH 08/11] shorten some sections --- cookbook/exploring_rbfe_settings.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index 80c307c..577e8b4 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -17,7 +17,7 @@ "source": [ "## Sampling and performance settings\n", "\n", - "This section covers the settings that control how much sampling is performed and therefore the computational cost and precision of a calculation.\n", + "This section covers the settings that affect sampling, computational cost and statistical precision.\n", "\n", "### Inspecting the default settings\n", "\n", @@ -410,7 +410,7 @@ "\n", "#### Changing the production length\n", "\n", - "The following example shows how to change the production length for a `RelativeHybridTopologyProtocol` object. In this example, we will double the production length from `5 ns` to `10 ns` per lambda window. Note the use of units from the `openff.units` package to specify the value is in nanoseconds. We will then calculate the new total simulation time and compare it to the default settings.\n" + "The following example shows how to double the production simulation length for the `RelativeHybridTopologyProtocol` and then calculates the new total simulation time and compares it to the default settings.\n" ], "id": "771ba6b9" }, @@ -594,7 +594,7 @@ "\n", "#### Changing the number of lambda windows\n", "\n", - "The following example shows how to change the number of lambda windows for a `RelativeHybridTopologyProtocol` object. In this example, we will double the number of lambda windows from `11` to `22`. We will then calculate the new total simulation time and compare it to the default settings. Note the `lambda_windows` setting only accepts integer values and generates a uniform spacing of lambda values between `0` and `1`.\n" + "The following example shows how to double the number of lambda windows for the `RelativeHybridTopologyProtocol` to `22` and calculates the new total simulation time. Note the `lambda_windows` setting only accepts integer values and generates a uniform spacing of lambda values between `0` and `1`.\n" ], "id": "3d6bb208a91ce6d0" }, From 5d9d93909663c6e1a4daf170d4c703e3ad8bc7c9 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Mon, 26 Oct 2026 13:19:23 +0000 Subject: [PATCH 09/11] more cleanup --- cookbook/exploring_rbfe_settings.ipynb | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index 577e8b4..5e41305 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -21,7 +21,7 @@ "\n", "### Inspecting the default settings\n", "\n", - "The default settings for each of the protocols are designed to be a good starting point for most alchemical transformations with the goal of obtaining reliable estimates with sufficient precision for decision-making at a reasonable computational cost. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", + "The default settings for each of the protocols are designed to be a starting point for most alchemical transformations with the goal of obtaining reliable estimates with sufficient precision for decision-making at a reasonable computational cost. However, depending on the specific system and transformation being studied, you may want to adjust these settings to improve convergence or precision of the estimation.\n", "\n", "The default settings for the `RelativeHybridTopologyProtocol` can be inspected by:" ] @@ -52,15 +52,14 @@ "\n", "\n", "\n", - "The main parameters that determine the amount of sampling performed and the associated cost for these protocols are:\n", + "The main parameters that determine the amount of sampling performed and the associated seeing for the `RelativeHybridTopologyProtocol` are:\n", "\n", - "| Setting | Description | `RelativeHybridTopologyProtocol` | `AbsoluteBindingProtocol` and `SepTopProtocol` |\n", - "|---------------------------------|--------------------------------------------------------------------------------------------------------------------------|--------------------------------------------|--------------------------------------------------------------|\n", - "| **Lambda windows** | The number of discrete alchemical states simulated during the transformation. | `lambda_settings.lambda_windows` | `{complex,solvent}_lambda_settings` |\n", - "| **Production length** | The length of the production simulation performed at each lambda window. | `simulation_settings.production_length` | `{complex,solvent}_simulation_settings.production_length` |\n", - "| **Equilibration length** | The length of the equilibration simulation performed at each lambda window using the selected multisate sampling method. | `simulation_settings.equilibration_length` | `{complex,solvent}_simulation_settings.equilibration_length` |\n", - "| **Pre-equilibration MD length** | The length of the single MD simulation per leg used to determine the restraint atoms in complex legs. | Not applicable | `{complex,solvent}_equil_simulation_settings` |\n", - "| **Protocol repeats** | The number of independent executions of the protocol. | `protocol_repeats` | `protocol_repeats` |\n", + "| Setting | Description | Setting |\n", + "|---------------------------------|--------------------------------------------------------------------------------------------------------------------------|--------------------------------------------|\n", + "| **Lambda windows** | The number of discrete alchemical states simulated during the transformation. | `lambda_settings.lambda_windows` |\n", + "| **Production length** | The length of the production simulation performed at each lambda window. | `simulation_settings.production_length` |\n", + "| **Equilibration length** | The length of the equilibration simulation performed at each lambda window using the selected multisate sampling method. | `simulation_settings.equilibration_length` | |\n", + "| **Protocol repeats** | The number of independent executions of the protocol. | `protocol_repeats` |\n", "\n", "The appropriate changes to these settings are transformation dependent with no single setting guaranteeing sufficient convergence or accuracy in general. These settings provide a good starting point to investigate the affect of additional sampling on the transformation under study however, but will also increase the computational cost of the simulation. The following example retrieves them from the default settings:" ], @@ -287,11 +286,11 @@ "source": [ "### Estimating the cost of the ABFE and SepTop protocols\n", "\n", - "The `AbsoluteBindingProtocol` and `SepTopProtocol` have a slightly different cost profile. Before the production simulation, each leg first runs a single, longer MD equilibration (NVT equilibration, NPT equilibration and a production run) that is used to determine the restraint atoms for the complex leg and seed the solvent leg. This simulation is run once per leg and repeat, **not** once per lambda window, so the total simulation time becomes:\n", + "The `AbsoluteBindingProtocol` and `SepTopProtocol` have a different cost profile. Each leg runs a single, longer MD equilibration (NVT equilibration, NPT equilibration and a production run) to determine the restraint atoms for the complex leg and seed the complex and solvent production simulations. This simulation is run once per leg and repeat, **not** once per lambda window, so the total simulation time becomes:\n", "\n", "$$T_{total} = N_{repeats} \\times \\sum_{legs} \\left[ T_{MD} + N_{\\lambda windows} \\times (T_{production} + T_{equilibration}) \\right]$$\n", "\n", - "where $T_{MD}$ is the sum of the NVT equilibration, equilibration and production lengths of the leg's `*_equil_simulation_settings`, and the remaining terms come from its `*_simulation_settings`. Both protocols require the number of replicas to match the number of lambda windows, so `n_replicas` gives the number of windows. The legs are `complex` and `solvent`." + "where $T_{MD}$ is the sum of the NVT equilibration, equilibration and production lengths of the leg's `*_equil_simulation_settings`, and the remaining terms come from its `*_simulation_settings`. The legs are `complex` and `solvent`." ], "id": "5017f2cf" }, From b669049f5bd3b7d312578270b1eb4b20a84ae66c Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Mon, 26 Oct 2026 16:28:58 +0000 Subject: [PATCH 10/11] refactor to group other protocols at the end --- cookbook/exploring_rbfe_settings.ipynb | 684 +++++++++++++------------ 1 file changed, 345 insertions(+), 339 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index 5e41305..f1b9018 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -103,370 +103,358 @@ "metadata": {}, "cell_type": "markdown", "source": [ - "#### Other protocols\n", + "### Estimating the computational cost of a protocol\n", + "\n", + "A useful first approximation of the computational cost of a protocol execution is the total simulation time in nanoseconds across all lambda windows and repeats. This can be calculated as:\n", "\n", - "Similar settings exist for the `AbsoluteBindingProtocol` and `SepTopProtocol`, but they are stored per leg (`complex_*` and `solvent_*`) as those protocols execute both `complex` and `solvent` legs within a single transformation. In addition, each leg has a `*_equil_simulation_settings` object controlling an initial MD simulation used to determine the restraint atoms in the complex legs, which is described in the cost estimate below.\n", + "$$T_{total} = N_{repeats} \\times N_{\\lambda windows} \\times (T_{production} + T_{equilibration})$$\n", "\n", - "**AbsoluteBindingProtocol**" + "So for the default settings the total simulation time is:" ], - "id": "a1b2c3d4e5f60001" + "id": "704d5dc6e6f895f5" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.502740Z", - "start_time": "2026-10-26T11:45:14.398886Z" + "end_time": "2026-10-26T11:45:16.681723Z", + "start_time": "2026-10-26T11:45:16.656346Z" } }, "cell_type": "code", "source": [ - "from openfe.protocols.openmm_afe import AbsoluteBindingProtocol\n", - "\n", - "abfe_default_settings = AbsoluteBindingProtocol.default_settings()\n", - "\n", - "print(f\"Independent repeats: {abfe_default_settings.protocol_repeats}\")\n", - "for leg in [\"complex\", \"solvent\"]:\n", - " equil = getattr(abfe_default_settings, f\"{leg}_equil_simulation_settings\")\n", - " sampling = getattr(abfe_default_settings, f\"{leg}_simulation_settings\")\n", + "def estimate_rbfe_cost(settings):\n", + " n_repeats = settings.protocol_repeats\n", + " n_windows = settings.lambda_settings.lambda_windows\n", + " production_length = settings.simulation_settings.production_length\n", + " equilibration_length = settings.simulation_settings.equilibration_length\n", "\n", - " print(f\"{leg}:\")\n", - " print(f\" Lambda windows: {sampling.n_replicas}\")\n", - " print(f\" Production length per state: {sampling.production_length.m} ns\")\n", - " print(f\" Equilibration length per state: {sampling.equilibration_length.m} ns\")\n", - " print(f\" Pre-equilibration MD length: {(equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length).m} ns\")" + " total_sim_time = n_repeats * n_windows * (production_length + equilibration_length)\n", + " return total_sim_time\n", + "total_sim_time_per_leg = estimate_rbfe_cost(settings)\n", + "print(f\"Total simulation time: {total_sim_time_per_leg.m} ns\")" ], - "id": "a1b2c3d4e5f60002", + "id": "fd81e13e2c34684d", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Independent repeats: 3\n", - "complex:\n", - " Lambda windows: 30\n", - " Production length per state: 10.0 ns\n", - " Equilibration length per state: 1 ns\n", - " Pre-equilibration MD length: 5.75 ns\n", - "solvent:\n", - " Lambda windows: 14\n", - " Production length per state: 10.0 ns\n", - " Equilibration length per state: 1.0 ns\n", - " Pre-equilibration MD length: 0.8 ns\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/Users/joshua/mambaforge/envs/openfe_dev/lib/python3.12/site-packages/Bio/Application/__init__.py:39: BiopythonDeprecationWarning: The Bio.Application modules and modules relying on it have been deprecated.\n", - "\n", - "Due to the on going maintenance burden of keeping command line application\n", - "wrappers up to date, we have decided to deprecate and eventually remove these\n", - "modules.\n", - "\n", - "We instead now recommend building your command line and invoking it directly\n", - "with the subprocess module.\n", - " warnings.warn(\n" + "Total simulation time: 198.0 ns\n" ] } ], - "execution_count": 3 + "execution_count": 5 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "**SepTopProtocol**" + "The `RelativeHybridTopologyProtocol` simulates the complex and solvent legs of a transformation independently, so if both use the same settings the total cost of the protocol is double the above estimate. It is useful to consider this estimate when deciding how to adjust the settings for a given transformation." + ], + "id": "3ab10c2ecebbd3df" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Increasing simulation length\n", + "\n", + "This section describes how to increase the amount of production simulation performed at each lambda window and estimate the increase in computational cost.\n", + "\n", + "
\n", + "Warning: Extending an existing simulation is not currently supported. Increasing the simulation length after running with default settings means the protocol must be run again from the start with the new settings, so the full cost of the longer simulation is incurred.\n", + "
\n", + "\n", + "#### When to increase simulation length\n", + "\n", + "Consider increasing the simulation length **before** running the simulation when:\n", + "\n", + "- The transformation involves large structural changes due to for example changes in formal charge. In this case you may also need to consider increasing the number of lambda windows see the next section for more information.\n", + "\n", + "Consider increasing the simulation length **after** running the simulation when:\n", + "\n", + "- Free-energy estimates differ substantially between independent repeats of the protocol.\n", + "- Forward and reverse free-energy estimates disagree by more than 1kbT at over 50% of the simulation length and their uncertainties do not overlap.[[1]](https://livecomsjournal.org/index.php/livecoms/article/view/v2i1e18378)\n", + "- A large number of samples are discarded during analysis due to the automated equilibration and decorrelation methods.\n", + "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows, you may also need to consider increasing the number of lambda windows see the next section for more information.\n", + "\n", + "This information is available through the automated analysis performed by the protocols, see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for the `RelativeHybridTopologyProtocol` for more information. It is recommended to review the analysis metrics before deciding to increase the simulation length.\n", + "\n", + "#### When not to increase simulation length\n", + "\n", + "Increasing the simulation length is not always the best solution. For example changes to the number of lambda windows or the spacing between them can sometimes be more effective when you have poor overlap between neighboring lambda windows and in some cases a combination of both maybe required.\n", + "\n", + "#### Changing the production length\n", + "\n", + "The following example shows how to double the production simulation length for the `RelativeHybridTopologyProtocol` and then calculates the new total simulation time and compares it to the default settings.\n" ], - "id": "a1b2c3d4e5f60003" + "id": "771ba6b9" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.632171Z", - "start_time": "2026-10-26T11:45:16.525149Z" + "end_time": "2026-10-26T11:45:16.910948Z", + "start_time": "2026-10-26T11:45:16.853625Z" } }, "cell_type": "code", "source": [ - "from openfe.protocols.openmm_septop import SepTopProtocol\n", + "from openff.units import unit\n", "\n", - "septop_default_settings = SepTopProtocol.default_settings()\n", + "long_settings = RelativeHybridTopologyProtocol.default_settings()\n", "\n", - "print(f\"Independent repeats: {septop_default_settings.protocol_repeats}\")\n", - "for leg in [\"complex\", \"solvent\"]:\n", - " equil = getattr(septop_default_settings, f\"{leg}_equil_simulation_settings\")\n", - " sampling = getattr(septop_default_settings, f\"{leg}_simulation_settings\")\n", + "# Increase the production length to 10 ns\n", + "long_settings.simulation_settings.production_length = 10 * unit.nanoseconds\n", "\n", - " print(f\"{leg}:\")\n", - " print(f\" Lambda windows: {sampling.n_replicas}\")\n", - " print(f\" Production length per state: {sampling.production_length.m} ns\")\n", - " print(f\" Equilibration length per state: {sampling.equilibration_length.m} ns\")\n", - " print(f\" Pre-equilibration MD length: {(equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length).m} ns\")" + "# Create a new protocol object with the updated settings\n", + "protocol = RelativeHybridTopologyProtocol(settings=long_settings)\n", + "\n", + "print(f\"Updated production length: {protocol.settings.simulation_settings.production_length.m} ns\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_total_sim_time_per_leg = estimate_rbfe_cost(protocol.settings)\n", + "print(f\"New total simulation time: {new_total_sim_time_per_leg.m} ns\")\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_total_sim_time_per_leg - total_sim_time_per_leg) / total_sim_time_per_leg) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" ], - "id": "a1b2c3d4e5f60004", + "id": "2cf9b3c7030df36f", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Independent repeats: 3\n", - "complex:\n", - " Lambda windows: 19\n", - " Production length per state: 10.0 ns\n", - " Equilibration length per state: 1.0 ns\n", - " Pre-equilibration MD length: 2.2 ns\n", - "solvent:\n", - " Lambda windows: 27\n", - " Production length per state: 10.0 ns\n", - " Equilibration length per state: 1.0 ns\n", - " Pre-equilibration MD length: 2.2 ns\n" + "Updated production length: 10 ns\n", + "New total simulation time: 363.0 ns\n", + "Percentage increase in total simulation time: 83.33%\n" ] } ], - "execution_count": 4 + "execution_count": 8 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "### Estimating the computational cost of a protocol\n", + "### Increasing the number of lambda windows\n", "\n", - "A useful first approximation of the computational cost of a protocol execution is the total simulation time in nanoseconds across all lambda windows and repeats. This can be calculated as:\n", + "Increasing the number of lambda windows introduces additional intermediate alchemical states along the transformation which can make the pathway easier to traverse and improve the overlap between neighboring states. This section describes how to increase the number of lambda windows and estimate the increase in computational cost.\n", "\n", - "$$T_{total} = N_{repeats} \\times N_{\\lambda windows} \\times (T_{production} + T_{equilibration})$$\n", + "#### When to increase the number of lambda windows\n", "\n", - "So for the default settings the total simulation time is:" + "Consider increasing the number of lambda windows **before** running the simulation when:\n", + "\n", + "- The transformation involves a change in net formal charge or substantial redistribution of charge, you may also need to consider increasing the simulation length see the previous section for more information.\n", + "\n", + "Consider increasing the number of lambda windows **after** running the simulation when:\n", + "\n", + "- Neighboring lambda windows have poor phase-space overlap or replica mixing statistics.\n", + "\n", + "#### When not to increase the number of lambda windows\n", + "\n", + "Increasing the number of lambda windows is not always the best solution, slow changes in ligand conformation or binding site rearrangements may require longer simulation lengths to achieve convergence.\n", + "\n", + "#### Changing the number of lambda windows\n", + "\n", + "The following example shows how to double the number of lambda windows for the `RelativeHybridTopologyProtocol` to `22` and calculates the new total simulation time. Note the `lambda_windows` setting only accepts integer values and generates a uniform spacing of lambda values between `0` and `1`.\n" ], - "id": "704d5dc6e6f895f5" + "id": "3d6bb208a91ce6d0" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.681723Z", - "start_time": "2026-10-26T11:45:16.656346Z" + "end_time": "2026-10-26T11:45:17.161411Z", + "start_time": "2026-10-26T11:45:17.043534Z" } }, "cell_type": "code", "source": [ - "def estimate_rbfe_cost(settings):\n", - " n_repeats = settings.protocol_repeats\n", - " n_windows = settings.lambda_settings.lambda_windows\n", - " production_length = settings.simulation_settings.production_length\n", - " equilibration_length = settings.simulation_settings.equilibration_length\n", + "more_window_settings = RelativeHybridTopologyProtocol.default_settings()\n", "\n", - " total_sim_time = n_repeats * n_windows * (production_length + equilibration_length)\n", - " return total_sim_time\n", - "total_sim_time_per_leg = estimate_rbfe_cost(settings)\n", - "print(f\"Total simulation time: {total_sim_time_per_leg.m} ns\")" + "# Increase the number of lambda windows to 22\n", + "more_window_settings.lambda_settings.lambda_windows = 22\n", + "# The number of replicas must be increased to match the number of lambda windows\n", + "more_window_settings.simulation_settings.n_replicas = 22\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "protocol = RelativeHybridTopologyProtocol(settings=more_window_settings)\n", + "\n", + "print(f\"Updated number of lambda windows: {protocol.settings.lambda_settings.lambda_windows}\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_total_sim_time_per_leg = estimate_rbfe_cost(protocol.settings)\n", + "print(f\"New total simulation time: {new_total_sim_time_per_leg.m} ns\")\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_total_sim_time_per_leg - total_sim_time_per_leg) / total_sim_time_per_leg) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" ], - "id": "fd81e13e2c34684d", + "id": "a8a0a2b5ba644eb4", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Total simulation time: 198.0 ns\n" + "Updated number of lambda windows: 22\n", + "New total simulation time: 396.0 ns\n", + "Percentage increase in total simulation time: 100.00%\n" ] } ], - "execution_count": 5 - }, - { - "metadata": {}, - "cell_type": "markdown", - "source": [ - "The `RelativeHybridTopologyProtocol` simulates the complex and solvent legs of a transformation independently, so if both use the same settings the total cost of the protocol is double the above estimate. It is useful to consider this estimate when deciding how to adjust the settings for a given transformation." - ], - "id": "3ab10c2ecebbd3df" + "execution_count": 11 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "### Estimating the cost of the ABFE and SepTop protocols\n", + "### Other protocols\n", "\n", - "The `AbsoluteBindingProtocol` and `SepTopProtocol` have a different cost profile. Each leg runs a single, longer MD equilibration (NVT equilibration, NPT equilibration and a production run) to determine the restraint atoms for the complex leg and seed the complex and solvent production simulations. This simulation is run once per leg and repeat, **not** once per lambda window, so the total simulation time becomes:\n", + "The above concepts apply to most free energy protocols, but the settings are located differently. Here we will briefly show how to inspect, and change them for the following protocols:\n", + " - `AbsoluteBindingProtocol`\n", + " - `SepTopProtocol`\n", "\n", - "$$T_{total} = N_{repeats} \\times \\sum_{legs} \\left[ T_{MD} + N_{\\lambda windows} \\times (T_{production} + T_{equilibration}) \\right]$$\n", - "\n", - "where $T_{MD}$ is the sum of the NVT equilibration, equilibration and production lengths of the leg's `*_equil_simulation_settings`, and the remaining terms come from its `*_simulation_settings`. The legs are `complex` and `solvent`." + "These protocols expose similar settings, but they are stored per leg (`complex_*` and `solvent_*`) as those protocols execute both `complex` and `solvent` legs within a single transformation. In addition, each leg has a `*_equil_simulation_settings` object controlling an initial MD simulation used to determine the restraint atoms in the complex legs and seed the solvent legs, which is described in the cost estimate for the `AbsoluteBindingProtocol`." ], - "id": "5017f2cf" + "id": "c1d2e3f4a5b60001" }, { "metadata": {}, "cell_type": "markdown", "source": [ - "**AbsoluteBindingProtocol**" + "### AbsoluteBindingProtocol\n", + "\n", + "#### Inspecting the default settings" ], - "id": "50aaaf76" + "id": "c1d2e3f4a5b60002" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.787042Z", - "start_time": "2026-10-26T11:45:16.711856Z" + "end_time": "2026-10-26T11:45:16.502740Z", + "start_time": "2026-10-26T11:45:14.398886Z" } }, "cell_type": "code", "source": [ - "from openff.units import unit\n", - "\n", - "def estimate_abfe_septop_cost(settings):\n", - " abfe_total_sim_time = 0 * unit.nanoseconds\n", - " for leg in [\"complex\", \"solvent\"]:\n", - " equil = getattr(settings, f\"{leg}_equil_simulation_settings\")\n", - " sampling = getattr(settings, f\"{leg}_simulation_settings\")\n", + "from openfe.protocols.openmm_afe import AbsoluteBindingProtocol\n", "\n", - " md_time = equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length\n", - " replica_exchange_time = sampling.n_replicas * (sampling.production_length + sampling.equilibration_length)\n", - " leg_time = md_time + replica_exchange_time\n", - " print(f\"{leg}: {md_time.m} ns MD + {replica_exchange_time.m} ns replica exchange = {leg_time.m} ns\")\n", - " abfe_total_sim_time += leg_time\n", + "abfe_default_settings = AbsoluteBindingProtocol.default_settings()\n", "\n", - " abfe_total_sim_time *= settings.protocol_repeats\n", - " print(f\"Total simulation time: {abfe_total_sim_time.m:.2f} ns\")\n", - " return abfe_total_sim_time\n", + "print(f\"Independent repeats: {abfe_default_settings.protocol_repeats}\")\n", + "for leg in [\"complex\", \"solvent\"]:\n", + " equil = getattr(abfe_default_settings, f\"{leg}_equil_simulation_settings\")\n", + " sampling = getattr(abfe_default_settings, f\"{leg}_simulation_settings\")\n", "\n", - "default_abfe_cost = estimate_abfe_septop_cost(abfe_default_settings)" + " print(f\"{leg}:\")\n", + " print(f\" Lambda windows: {sampling.n_replicas}\")\n", + " print(f\" Production length per state: {sampling.production_length.m} ns\")\n", + " print(f\" Equilibration length per state: {sampling.equilibration_length.m} ns\")\n", + " print(f\" Pre-equilibration MD length: {(equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length).m} ns\")" ], - "id": "0f8fb861", + "id": "a1b2c3d4e5f60002", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "complex: 5.75 ns MD + 330.0 ns replica exchange = 335.75 ns\n", - "solvent: 0.8 ns MD + 154.0 ns replica exchange = 154.8 ns\n", - "Total simulation time: 1471.65 ns\n" + "Independent repeats: 3\n", + "complex:\n", + " Lambda windows: 30\n", + " Production length per state: 10.0 ns\n", + " Equilibration length per state: 1 ns\n", + " Pre-equilibration MD length: 5.75 ns\n", + "solvent:\n", + " Lambda windows: 14\n", + " Production length per state: 10.0 ns\n", + " Equilibration length per state: 1.0 ns\n", + " Pre-equilibration MD length: 0.8 ns\n" ] - } - ], - "execution_count": 6 - }, - { - "metadata": {}, - "cell_type": "markdown", - "source": [ - "**SepTopProtocol**" - ], - "id": "b004911d" - }, - { - "metadata": { - "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.850622Z", - "start_time": "2026-10-26T11:45:16.790688Z" - } - }, - "cell_type": "code", - "source": "default_septop_cost = estimate_abfe_septop_cost(septop_default_settings)", - "id": "85ed4076", - "outputs": [ + }, { - "name": "stdout", + "name": "stderr", "output_type": "stream", "text": [ - "complex: 2.2 ns MD + 209.0 ns replica exchange = 211.2 ns\n", - "solvent: 2.2 ns MD + 297.0 ns replica exchange = 299.2 ns\n", - "Total simulation time: 1531.20 ns\n" + "/Users/joshua/mambaforge/envs/openfe_dev/lib/python3.12/site-packages/Bio/Application/__init__.py:39: BiopythonDeprecationWarning: The Bio.Application modules and modules relying on it have been deprecated.\n", + "\n", + "Due to the on going maintenance burden of keeping command line application\n", + "wrappers up to date, we have decided to deprecate and eventually remove these\n", + "modules.\n", + "\n", + "We instead now recommend building your command line and invoking it directly\n", + "with the subprocess module.\n", + " warnings.warn(\n" ] } ], - "execution_count": 7 + "execution_count": 3 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "### Increasing simulation length\n", - "\n", - "This section describes how to increase the amount of production simulation performed at each lambda window and estimate the increase in computational cost.\n", + "#### Estimating the computational cost\n", "\n", - "
\n", - "Warning: Extending an existing simulation is not currently supported. Increasing the simulation length after running with default settings means the protocol must be run again from the start with the new settings, so the full cost of the longer simulation is incurred.\n", - "
\n", - "\n", - "#### When to increase simulation length\n", - "\n", - "Consider increasing the simulation length **before** running the simulation when:\n", + "The `AbsoluteBindingProtocol` and `SepTopProtocol` have a different cost profile to the `RelativeHybridTopologyProtocol`. Each leg runs a single, longer MD equilibration (NVT equilibration, NPT equilibration and a production run). This simulation is run once per leg and repeat, **not** once per lambda window, so the total simulation time becomes:\n", "\n", - "- The transformation involves large structural changes due to for example changes in formal charge. In this case you may also need to consider increasing the number of lambda windows see the next section for more information.\n", - "\n", - "Consider increasing the simulation length **after** running the simulation when:\n", - "\n", - "- Free-energy estimates differ substantially between independent repeats of the protocol.\n", - "- Forward and reverse free-energy estimates disagree by more than 1kbT at over 50% of the simulation length and their uncertainties do not overlap.[[1]](https://livecomsjournal.org/index.php/livecoms/article/view/v2i1e18378)\n", - "- A large number of samples are discarded during analysis due to the automated equilibration and decorrelation methods.\n", - "- Replica exchange or MBAR analysis indicates that alchemical overlap or replica mixing is poor between lambda windows, you may also need to consider increasing the number of lambda windows see the next section for more information.\n", - "\n", - "This information is available through the automated analysis performed by the protocols, see the list of available [analysis metrics](https://docs.openfree.energy/en/latest/guide/protocols/relativehybridtopology.html#analysis) for the `RelativeHybridTopologyProtocol` for more information. It is recommended to review the analysis metrics before deciding to increase the simulation length.\n", - "\n", - "#### When not to increase simulation length\n", - "\n", - "Increasing the simulation length is not always the best solution. For example changes to the number of lambda windows or the spacing between them can sometimes be more effective when you have poor overlap between neighboring lambda windows and in some cases a combination of both maybe required.\n", - "\n", - "#### Changing the production length\n", + "$$T_{total} = N_{repeats} \\times \\sum_{legs} \\left[ T_{MD} + N_{\\lambda windows} \\times (T_{production} + T_{equilibration}) \\right]$$\n", "\n", - "The following example shows how to double the production simulation length for the `RelativeHybridTopologyProtocol` and then calculates the new total simulation time and compares it to the default settings.\n" + "where $T_{MD}$ is the sum of the NVT equilibration, equilibration and production lengths of the leg's `*_equil_simulation_settings`, and the remaining terms come from its `*_simulation_settings`." ], - "id": "771ba6b9" + "id": "c1d2e3f4a5b60003" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.910948Z", - "start_time": "2026-10-26T11:45:16.853625Z" + "end_time": "2026-10-26T16:26:32.709914Z", + "start_time": "2026-10-26T16:26:32.426034Z" } }, "cell_type": "code", "source": [ - "from openff.units import unit\n", - "\n", - "long_settings = RelativeHybridTopologyProtocol.default_settings()\n", - "\n", - "# Increase the production length to 10 ns\n", - "long_settings.simulation_settings.production_length = 10 * unit.nanoseconds\n", - "\n", - "# Create a new protocol object with the updated settings\n", - "protocol = RelativeHybridTopologyProtocol(settings=long_settings)\n", + "def estimate_abfe_septop_cost(settings):\n", + " abfe_total_sim_time = 0 * unit.nanoseconds\n", + " for leg in [\"complex\", \"solvent\"]:\n", + " equil = getattr(settings, f\"{leg}_equil_simulation_settings\")\n", + " sampling = getattr(settings, f\"{leg}_simulation_settings\")\n", "\n", - "print(f\"Updated production length: {protocol.settings.simulation_settings.production_length.m} ns\")\n", + " md_time = equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length\n", + " replica_exchange_time = sampling.n_replicas * (sampling.production_length + sampling.equilibration_length)\n", + " leg_time = md_time + replica_exchange_time\n", + " print(f\"{leg}: {md_time.m} ns MD + {replica_exchange_time.m} ns replica exchange = {leg_time.m} ns\")\n", + " abfe_total_sim_time += leg_time\n", "\n", - "# Calculate the new total simulation time\n", - "new_total_sim_time_per_leg = estimate_rbfe_cost(protocol.settings)\n", - "print(f\"New total simulation time: {new_total_sim_time_per_leg.m} ns\")\n", + " abfe_total_sim_time *= settings.protocol_repeats\n", + " print(f\"Total simulation time: {abfe_total_sim_time.m:.2f} ns\")\n", + " return abfe_total_sim_time\n", "\n", - "# Print the percentage increase in total simulation time\n", - "percentage_increase = ((new_total_sim_time_per_leg - total_sim_time_per_leg) / total_sim_time_per_leg) * 100\n", - "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + "default_abfe_cost = estimate_abfe_septop_cost(abfe_default_settings)" ], - "id": "2cf9b3c7030df36f", + "id": "0f8fb861", "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Updated production length: 10 ns\n", - "New total simulation time: 363.0 ns\n", - "Percentage increase in total simulation time: 83.33%\n" + "ename": "NameError", + "evalue": "name 'abfe_default_settings' is not defined", + "output_type": "error", + "traceback": [ + "\u001B[31m---------------------------------------------------------------------------\u001B[39m", + "\u001B[31mNameError\u001B[39m Traceback (most recent call last)", + "\u001B[36mCell\u001B[39m\u001B[36m \u001B[39m\u001B[32mIn[1]\u001B[39m\u001B[32m, line 17\u001B[39m\n\u001B[32m 14\u001B[39m \u001B[38;5;28mprint\u001B[39m(\u001B[33mf\u001B[39m\u001B[33m\"\u001B[39m\u001B[33mTotal simulation time: \u001B[39m\u001B[38;5;132;01m{\u001B[39;00mabfe_total_sim_time.m\u001B[38;5;132;01m:\u001B[39;00m\u001B[33m.2f\u001B[39m\u001B[38;5;132;01m}\u001B[39;00m\u001B[33m ns\u001B[39m\u001B[33m\"\u001B[39m)\n\u001B[32m 15\u001B[39m \u001B[38;5;28;01mreturn\u001B[39;00m abfe_total_sim_time\n\u001B[32m---> \u001B[39m\u001B[32m17\u001B[39m default_abfe_cost = estimate_abfe_septop_cost(\u001B[43mabfe_default_settings\u001B[49m)\n", + "\u001B[31mNameError\u001B[39m: name 'abfe_default_settings' is not defined" ] } ], - "execution_count": 8 + "execution_count": 1 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "##### Other protocols\n", - "\n", - "The same idea applies to each of the other protocols, only the location of the setting changes. Protocols that simulate several legs (e.g. complex and solvent) store the settings per leg, so you must change each leg to be extended individually.\n", + "#### Increasing simulation length\n", "\n", - "**AbsoluteBindingProtocol**" + "The production length is set per leg, so each leg to be extended must be changed individually." ], - "id": "d4f1a7c2e9b34a60" + "id": "c1d2e3f4a5b60004" }, { "metadata": { @@ -519,225 +507,243 @@ "metadata": {}, "cell_type": "markdown", "source": [ - "**SepTopProtocol**" + "#### Changing the number of lambda windows\n", + "\n", + "The `AbsoluteBindingProtocol` and `SepTopProtocol` both allow for exact control of the lambda schedule of each leg requiring an explicitly defined list of values, one entry per window. All of the lists for a leg must have the same length, and `n_replicas` in the matching `*_simulation_settings` must be set to the number of windows.\n", + "\n", + "The schedule is defined by `lambda_elec`, `lambda_vdw` and `lambda_restraints`." ], - "id": "f2b3c4d5e6a71829" + "id": "c1d2e3f4a5b60005" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.016519Z", - "start_time": "2026-10-26T11:45:16.983332Z" + "end_time": "2026-10-26T11:45:17.218775Z", + "start_time": "2026-10-26T11:45:17.183233Z" } }, "cell_type": "code", "source": [ - "from openfe.protocols.openmm_septop import SepTopProtocol\n", + "abfe_window_settings = AbsoluteBindingProtocol.default_settings()\n", "\n", - "septop_long_settings = SepTopProtocol.default_settings()\n", + "# Define the new lambda schedule for the complex leg by adding 3 more charge windows\n", + "abfe_window_settings.complex_lambda_settings.lambda_elec = [0.0] * 6 + [0.1, 0.2, 0.3, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.85, 0.9] + [1.0] * 15\n", + "# Extra vdW windows are added to match the length of the new charge schedule\n", + "abfe_window_settings.complex_lambda_settings.lambda_vdw = [0.0] * 19 + [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]\n", + "# Extra restraint windows are added to match the length of the new charge schedule\n", + "abfe_window_settings.complex_lambda_settings.lambda_restraints = [0.0, 0.2, 0.4, 0.6, 0.8] + [1.0] * 28\n", "\n", - "septop_long_settings.complex_simulation_settings.production_length = 20 * unit.nanoseconds\n", - "septop_long_settings.solvent_simulation_settings.production_length = 20 * unit.nanoseconds\n", + "# The number of replicas must match the number of lambda windows\n", + "abfe_window_settings.complex_simulation_settings.n_replicas = len(abfe_window_settings.complex_lambda_settings.lambda_vdw)\n", "\n", "# Create a new protocol object with the updated settings\n", - "septop_protocol = SepTopProtocol(settings=septop_long_settings)\n", - "\n", - "print(f\"Updated complex production length: {septop_protocol.settings.complex_simulation_settings.production_length.m} ns\")\n", - "print(f\"Updated solvent production length: {septop_protocol.settings.solvent_simulation_settings.production_length.m} ns\")\n", + "abfe_protocol = AbsoluteBindingProtocol(settings=abfe_window_settings)\n", "\n", "# Calculate the new total simulation time\n", - "new_septop_cost = estimate_abfe_septop_cost(septop_long_settings)\n", + "new_abfe_cost = estimate_abfe_septop_cost(abfe_window_settings)\n", "\n", "# Print the percentage increase in total simulation time\n", - "percentage_increase = ((new_septop_cost - default_septop_cost) / default_septop_cost) * 100\n", + "percentage_increase = ((new_abfe_cost - default_abfe_cost) / default_abfe_cost) * 100\n", "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" ], - "id": "a3c4d5e6f7b82930", + "id": "b1c2d3e4f5a60002", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Updated complex production length: 20 ns\n", - "Updated solvent production length: 20 ns\n", - "complex: 2.2 ns MD + 399.0 ns replica exchange = 401.2 ns\n", - "solvent: 2.2 ns MD + 567.0 ns replica exchange = 569.2 ns\n", - "Total simulation time: 2911.20 ns\n", - "Percentage increase in total simulation time: 90.13%\n" + "complex: 5.75 ns MD + 363.0 ns replica exchange = 368.75 ns\n", + "solvent: 0.8 ns MD + 154.0 ns replica exchange = 154.8 ns\n", + "Total simulation time: 1570.65 ns\n", + "Percentage increase in total simulation time: 6.73%\n" ] } ], - "execution_count": 10 + "execution_count": 12 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "### Increasing the number of lambda windows\n", - "\n", - "Increasing the number of lambda windows introduces additional intermediate alchemical states along the transformation which can make the pathway easier to traverse and improve the overlap between neighboring states. This section describes how to increase the number of lambda windows and estimate the increase in computational cost.\n", - "\n", - "#### When to increase the number of lambda windows\n", - "\n", - "Consider increasing the number of lambda windows **before** running the simulation when:\n", - "\n", - "- The transformation involves a change in net formal charge or substantial redistribution of charge, you may also need to consider increasing the simulation length see the previous section for more information.\n", - "\n", - "Consider increasing the number of lambda windows **after** running the simulation when:\n", - "\n", - "- Neighboring lambda windows have poor phase-space overlap or replica mixing statistics.\n", - "\n", - "#### When not to increase the number of lambda windows\n", - "\n", - "Increasing the number of lambda windows is not always the best solution, slow changes in ligand conformation or binding site rearrangements may require longer simulation lengths to achieve convergence.\n", - "\n", - "#### Changing the number of lambda windows\n", - "\n", - "The following example shows how to double the number of lambda windows for the `RelativeHybridTopologyProtocol` to `22` and calculates the new total simulation time. Note the `lambda_windows` setting only accepts integer values and generates a uniform spacing of lambda values between `0` and `1`.\n" + "
\n", + "Tip: Changing the number of lambda windows or the lambda schedule for the AbsoluteBindingProtocol and SepTopProtocol by hand can be error prone. The protocols provide a validation method to check that the lambda schedules are valid for the leg in question, use them to validate your changes and catch issues before execution. See the example below.\n", + "
" ], - "id": "3d6bb208a91ce6d0" + "id": "2290219ca360bbd3" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.161411Z", - "start_time": "2026-10-26T11:45:17.043534Z" + "end_time": "2026-10-26T11:45:17.251783Z", + "start_time": "2026-10-26T11:45:17.234232Z" } }, "cell_type": "code", "source": [ - "more_window_settings = RelativeHybridTopologyProtocol.default_settings()\n", - "\n", - "# Increase the number of lambda windows to 22\n", - "more_window_settings.lambda_settings.lambda_windows = 22\n", - "# The number of replicas must be increased to match the number of lambda windows\n", - "more_window_settings.simulation_settings.n_replicas = 22\n", + "# Validate the lambda schedules for the complex leg\n", + "AbsoluteBindingProtocol._validate_lambda_schedule(\n", + " abfe_window_settings.complex_lambda_settings,\n", + " abfe_window_settings.complex_simulation_settings\n", + ")" + ], + "id": "6c9de9ae4d246d9c", + "outputs": [], + "execution_count": 13 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### SepTopProtocol\n", "\n", - "# Create a new protocol object with the updated settings\n", - "protocol = RelativeHybridTopologyProtocol(settings=more_window_settings)\n", + "#### Inspecting the default settings" + ], + "id": "c1d2e3f4a5b60006" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T11:45:16.632171Z", + "start_time": "2026-10-26T11:45:16.525149Z" + } + }, + "cell_type": "code", + "source": [ + "from openfe.protocols.openmm_septop import SepTopProtocol\n", "\n", - "print(f\"Updated number of lambda windows: {protocol.settings.lambda_settings.lambda_windows}\")\n", + "septop_default_settings = SepTopProtocol.default_settings()\n", "\n", - "# Calculate the new total simulation time\n", - "new_total_sim_time_per_leg = estimate_rbfe_cost(protocol.settings)\n", - "print(f\"New total simulation time: {new_total_sim_time_per_leg.m} ns\")\n", + "print(f\"Independent repeats: {septop_default_settings.protocol_repeats}\")\n", + "for leg in [\"complex\", \"solvent\"]:\n", + " equil = getattr(septop_default_settings, f\"{leg}_equil_simulation_settings\")\n", + " sampling = getattr(septop_default_settings, f\"{leg}_simulation_settings\")\n", "\n", - "# Print the percentage increase in total simulation time\n", - "percentage_increase = ((new_total_sim_time_per_leg - total_sim_time_per_leg) / total_sim_time_per_leg) * 100\n", - "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + " print(f\"{leg}:\")\n", + " print(f\" Lambda windows: {sampling.n_replicas}\")\n", + " print(f\" Production length per state: {sampling.production_length.m} ns\")\n", + " print(f\" Equilibration length per state: {sampling.equilibration_length.m} ns\")\n", + " print(f\" Pre-equilibration MD length: {(equil.equilibration_length_nvt + equil.equilibration_length + equil.production_length).m} ns\")" ], - "id": "a8a0a2b5ba644eb4", + "id": "a1b2c3d4e5f60004", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Updated number of lambda windows: 22\n", - "New total simulation time: 396.0 ns\n", - "Percentage increase in total simulation time: 100.00%\n" + "Independent repeats: 3\n", + "complex:\n", + " Lambda windows: 19\n", + " Production length per state: 10.0 ns\n", + " Equilibration length per state: 1.0 ns\n", + " Pre-equilibration MD length: 2.2 ns\n", + "solvent:\n", + " Lambda windows: 27\n", + " Production length per state: 10.0 ns\n", + " Equilibration length per state: 1.0 ns\n", + " Pre-equilibration MD length: 2.2 ns\n" ] } ], - "execution_count": 11 + "execution_count": 4 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "##### Other protocols\n", - "\n", - "The `AbsoluteBindingProtocol` and `SepTopProtocol` allow for exact control of the lambda schedule of each leg requiring an explicitly defined list of values, one entry per window. All of the lists for a leg must have the same length, and `n_replicas` in the matching `*_simulation_settings` must be set to the number of windows. The examples below change the `complex` leg, the `solvent` leg can be changed in the same way using `solvent_lambda_settings` and `solvent_simulation_settings`.\n", + "#### Estimating the computational cost\n", "\n", - "**AbsoluteBindingProtocol**\n", - "\n", - "The schedule is defined by `lambda_elec`, `lambda_vdw` and `lambda_restraints`." + "The `SepTopProtocol` has the same cost profile as the `AbsoluteBindingProtocol`, so the `estimate_abfe_septop_cost` function defined above can be reused." ], - "id": "b1c2d3e4f5a60001" + "id": "c1d2e3f4a5b60007" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.218775Z", - "start_time": "2026-10-26T11:45:17.183233Z" + "end_time": "2026-10-26T11:45:16.850622Z", + "start_time": "2026-10-26T11:45:16.790688Z" } }, "cell_type": "code", - "source": [ - "abfe_window_settings = AbsoluteBindingProtocol.default_settings()\n", - "\n", - "# Define the new lambda schedule for the complex leg by adding 3 more charge windows\n", - "abfe_window_settings.complex_lambda_settings.lambda_elec = [0.0] * 6 + [0.1, 0.2, 0.3, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.85, 0.9] + [1.0] * 15\n", - "# Extra vdW windows are added to match the length of the new charge schedule\n", - "abfe_window_settings.complex_lambda_settings.lambda_vdw = [0.0] * 19 + [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]\n", - "# Extra restraint windows are added to match the length of the new charge schedule\n", - "abfe_window_settings.complex_lambda_settings.lambda_restraints = [0.0, 0.2, 0.4, 0.6, 0.8] + [1.0] * 28\n", - "\n", - "# The number of replicas must match the number of lambda windows\n", - "abfe_window_settings.complex_simulation_settings.n_replicas = len(abfe_window_settings.complex_lambda_settings.lambda_vdw)\n", - "\n", - "# Create a new protocol object with the updated settings\n", - "abfe_protocol = AbsoluteBindingProtocol(settings=abfe_window_settings)\n", - "\n", - "# Calculate the new total simulation time\n", - "new_abfe_cost = estimate_abfe_septop_cost(abfe_window_settings)\n", - "\n", - "# Print the percentage increase in total simulation time\n", - "percentage_increase = ((new_abfe_cost - default_abfe_cost) / default_abfe_cost) * 100\n", - "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" - ], - "id": "b1c2d3e4f5a60002", + "source": "default_septop_cost = estimate_abfe_septop_cost(septop_default_settings)", + "id": "85ed4076", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "complex: 5.75 ns MD + 363.0 ns replica exchange = 368.75 ns\n", - "solvent: 0.8 ns MD + 154.0 ns replica exchange = 154.8 ns\n", - "Total simulation time: 1570.65 ns\n", - "Percentage increase in total simulation time: 6.73%\n" + "complex: 2.2 ns MD + 209.0 ns replica exchange = 211.2 ns\n", + "solvent: 2.2 ns MD + 297.0 ns replica exchange = 299.2 ns\n", + "Total simulation time: 1531.20 ns\n" ] } ], - "execution_count": 12 + "execution_count": 7 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "
\n", - "Tip: Changing the number of lambda windows or the lambda schedule for the AbsoluteBindingProtocol and SepTopProtocol by hand can be error prone. The protocols provide a validation method to check that the lambda schedules are valid for the leg in question, use them to validate your changes and catch issues before execution. See the example below.\n", - "
" + "#### Increasing simulation length\n", + "\n", + "The production length is changed per leg in the same way as for the `AbsoluteBindingProtocol`." ], - "id": "2290219ca360bbd3" + "id": "c1d2e3f4a5b60008" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.251783Z", - "start_time": "2026-10-26T11:45:17.234232Z" + "end_time": "2026-10-26T11:45:17.016519Z", + "start_time": "2026-10-26T11:45:16.983332Z" } }, "cell_type": "code", "source": [ - "# Validate the lambda schedules for the complex leg\n", - "AbsoluteBindingProtocol._validate_lambda_schedule(\n", - " abfe_window_settings.complex_lambda_settings,\n", - " abfe_window_settings.complex_simulation_settings\n", - ")" + "from openfe.protocols.openmm_septop import SepTopProtocol\n", + "\n", + "septop_long_settings = SepTopProtocol.default_settings()\n", + "\n", + "septop_long_settings.complex_simulation_settings.production_length = 20 * unit.nanoseconds\n", + "septop_long_settings.solvent_simulation_settings.production_length = 20 * unit.nanoseconds\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "septop_protocol = SepTopProtocol(settings=septop_long_settings)\n", + "\n", + "print(f\"Updated complex production length: {septop_protocol.settings.complex_simulation_settings.production_length.m} ns\")\n", + "print(f\"Updated solvent production length: {septop_protocol.settings.solvent_simulation_settings.production_length.m} ns\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_septop_cost = estimate_abfe_septop_cost(septop_long_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_septop_cost - default_septop_cost) / default_septop_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" ], - "id": "6c9de9ae4d246d9c", - "outputs": [], - "execution_count": 13 + "id": "a3c4d5e6f7b82930", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated complex production length: 20 ns\n", + "Updated solvent production length: 20 ns\n", + "complex: 2.2 ns MD + 399.0 ns replica exchange = 401.2 ns\n", + "solvent: 2.2 ns MD + 567.0 ns replica exchange = 569.2 ns\n", + "Total simulation time: 2911.20 ns\n", + "Percentage increase in total simulation time: 90.13%\n" + ] + } + ], + "execution_count": 10 }, { "metadata": {}, "cell_type": "markdown", "source": [ - "**SepTopProtocol**\n", + "#### Changing the number of lambda windows\n", "\n", "The schedule is defined by the `A` and `B` states of each component: `lambda_elec_A`, `lambda_elec_B`, `lambda_vdw_A`, `lambda_vdw_B`, `lambda_restraints_A` and `lambda_restraints_B`. Like the `AbsoluteBindingProtocol` all the lambda value lists for a leg must be the same length." ], - "id": "b1c2d3e4f5a60003" + "id": "c1d2e3f4a5b60009" }, { "metadata": { From c1d1cd9c3ec2f4bf4571bd33c262e8c9e151b750 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Tue, 27 Oct 2026 09:23:47 +0000 Subject: [PATCH 11/11] add repeats and equilibration length updates --- cookbook/exploring_rbfe_settings.ipynb | 450 ++++++++++++++++++++++--- 1 file changed, 397 insertions(+), 53 deletions(-) diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb index f1b9018..f8dfd72 100644 --- a/cookbook/exploring_rbfe_settings.ipynb +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -29,8 +29,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:14.335951Z", - "start_time": "2026-10-26T11:45:01.020777Z" + "end_time": "2026-10-26T16:54:59.569668Z", + "start_time": "2026-10-26T16:54:45.868611Z" } }, "cell_type": "code", @@ -68,8 +68,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:14.395446Z", - "start_time": "2026-10-26T11:45:14.352357Z" + "end_time": "2026-10-26T16:54:59.598767Z", + "start_time": "2026-10-26T16:54:59.585534Z" } }, "cell_type": "code", @@ -116,8 +116,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.681723Z", - "start_time": "2026-10-26T11:45:16.656346Z" + "end_time": "2026-10-26T16:54:59.633570Z", + "start_time": "2026-10-26T16:54:59.618645Z" } }, "cell_type": "code", @@ -143,7 +143,7 @@ ] } ], - "execution_count": 5 + "execution_count": 3 }, { "metadata": {}, @@ -193,8 +193,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.910948Z", - "start_time": "2026-10-26T11:45:16.853625Z" + "end_time": "2026-10-26T16:54:59.664664Z", + "start_time": "2026-10-26T16:54:59.647568Z" } }, "cell_type": "code", @@ -231,7 +231,7 @@ ] } ], - "execution_count": 8 + "execution_count": 4 }, { "metadata": {}, @@ -264,8 +264,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.161411Z", - "start_time": "2026-10-26T11:45:17.043534Z" + "end_time": "2026-10-26T16:54:59.695800Z", + "start_time": "2026-10-26T16:54:59.679630Z" } }, "cell_type": "code", @@ -302,7 +302,139 @@ ] } ], - "execution_count": 11 + "execution_count": 5 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Increasing the number of independent repeats\n", + "\n", + "Independent repeats run the same protocol from the same input structures but different initial velocities. This is the default method used to quantify the uncertainty of a protocol estimate in **openfe** and is useful to assess the reproducibility of the result.\n", + "\n", + "#### When to run additional repeats\n", + "\n", + "Consider running additional repeats when:\n", + "\n", + "- You want to quantify the uncertainty of the free energy estimate beyond the estimator uncertainty provided by MBAR bootstrap resampling.\n", + "\n", + "#### When not to run additional repeats\n", + "\n", + "If the uncertainty of the free energy estimate is already sufficiently small for your application, additional repeats may not be necessary. The default number of repeats is designed to provide a reasonable estimate of uncertainty for most transformations.\n", + "\n", + "If the uncertainty is very large, it may be more efficient to investigate the cause of the variation between repeats which might be better effected by increasing the simulation length or number of lambda windows.\n", + "\n", + "#### Changing the number of repeats\n", + "\n", + "The following example shows how to increase the number of independent repeats per single execution of the `RelativeHybridTopologyProtocol` to `4`. Note that you can also execute the protocol multiple times with a single repeat per execution which may be more efficient, see the [quickrun guide](https://docs.openfree.energy/en/stable/guide/execution/quickrun_execution.html#parallel-execution-of-repeats-with-quickrun) for more information." + ], + "id": "26c868fb24073c7e" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T16:56:58.685724Z", + "start_time": "2026-10-26T16:56:58.647979Z" + } + }, + "cell_type": "code", + "source": [ + "more_repeats_settings = RelativeHybridTopologyProtocol.default_settings()\n", + "\n", + "# Increase the number of independent repeats to 4\n", + "more_repeats_settings.protocol_repeats = 4\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "protocol = RelativeHybridTopologyProtocol(settings=more_repeats_settings)\n", + "\n", + "print(f\"Updated number of independent repeats: {protocol.settings.protocol_repeats}\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_total_sim_time_per_leg = estimate_rbfe_cost(protocol.settings)\n", + "print(f\"New total simulation time: {new_total_sim_time_per_leg.m} ns\")\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_total_sim_time_per_leg - total_sim_time_per_leg) / total_sim_time_per_leg) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "8a8f3e35757daeb3", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated number of independent repeats: 4\n", + "New total simulation time: 264.0 ns\n", + "Percentage increase in total simulation time: 33.33%\n" + ] + } + ], + "execution_count": 19 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Increasing the equilibration length\n", + "\n", + "The equilibration simulation uses the selected multisate sampling method to relax the system at each lambda window before production sampling begins. No samples are collected during this phase for the final free energy estimate.\n", + "\n", + "#### When to increase the equilibration length\n", + "\n", + "Consider increasing the equilibration length **before** running the simulation when:\n", + "\n", + "- The input structures have large steric clashes or are far from equilibrium and the system is not expected to equilibrate quickly.\n", + "\n", + "#### When not to increase the equilibration length\n", + "\n", + "If the system is expected to require a long equilibration time, such those with membrane bound proteins, it may be more efficient to perform a longer equilibration simulation outside the protocol.\n", + "\n", + "#### Changing the equilibration length\n", + "\n", + "The following example shows how to double the equilibration simulation length for the `RelativeHybridTopologyProtocol` and calculates the new total simulation time. Note that the equilibration length is set per lambda window, so increasing it will increase the total simulation time linearly with the number of windows see increasing production length for more information on the associated cost." + ], + "id": "1c95943803da691e" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-26T16:56:13.313640Z", + "start_time": "2026-10-26T16:56:13.242224Z" + } + }, + "cell_type": "code", + "source": [ + "more_equilibration_settings = RelativeHybridTopologyProtocol.default_settings()\n", + "\n", + "# Increase the equilibration length to 2 ns\n", + "more_equilibration_settings.simulation_settings.equilibration_length = 2 * unit.nanoseconds\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "protocol = RelativeHybridTopologyProtocol(settings=more_equilibration_settings)\n", + "\n", + "print(f\"Updated equilibration length: {protocol.settings.simulation_settings.equilibration_length.m} ns\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_total_sim_time_per_leg = estimate_rbfe_cost(protocol.settings)\n", + "print(f\"New total simulation time: {new_total_sim_time_per_leg.m} ns\")\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_total_sim_time_per_leg - total_sim_time_per_leg) / total_sim_time_per_leg) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "5a8fdf02930d79e3", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated equilibration length: 2 ns\n", + "New total simulation time: 231.0 ns\n", + "Percentage increase in total simulation time: 16.67%\n" + ] + } + ], + "execution_count": 18 }, { "metadata": {}, @@ -314,7 +446,7 @@ " - `AbsoluteBindingProtocol`\n", " - `SepTopProtocol`\n", "\n", - "These protocols expose similar settings, but they are stored per leg (`complex_*` and `solvent_*`) as those protocols execute both `complex` and `solvent` legs within a single transformation. In addition, each leg has a `*_equil_simulation_settings` object controlling an initial MD simulation used to determine the restraint atoms in the complex legs and seed the solvent legs, which is described in the cost estimate for the `AbsoluteBindingProtocol`." + "As these protocols execute both the `complex` and `solvent` legs within a single transformation the settings are exposed with `complex_*` or `solvent_*` prefix. In addition, each leg has a `*_equil_simulation_settings` object controlling an initial MD simulation used to determine the restraint atoms in the complex legs and seed the solvent legs, which is described in the cost estimate for the `AbsoluteBindingProtocol`." ], "id": "c1d2e3f4a5b60001" }, @@ -331,8 +463,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.502740Z", - "start_time": "2026-10-26T11:45:14.398886Z" + "end_time": "2026-10-26T16:55:01.234321Z", + "start_time": "2026-10-26T16:54:59.797756Z" } }, "cell_type": "code", @@ -387,7 +519,7 @@ ] } ], - "execution_count": 3 + "execution_count": 8 }, { "metadata": {}, @@ -399,15 +531,15 @@ "\n", "$$T_{total} = N_{repeats} \\times \\sum_{legs} \\left[ T_{MD} + N_{\\lambda windows} \\times (T_{production} + T_{equilibration}) \\right]$$\n", "\n", - "where $T_{MD}$ is the sum of the NVT equilibration, equilibration and production lengths of the leg's `*_equil_simulation_settings`, and the remaining terms come from its `*_simulation_settings`." + "where $T_{MD}$ is the sum of the NVT equilibration. The equilibration and production lengths of the leg's are set by `*_equil_simulation_settings`, and the remaining terms come from its `*_simulation_settings`." ], "id": "c1d2e3f4a5b60003" }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T16:26:32.709914Z", - "start_time": "2026-10-26T16:26:32.426034Z" + "end_time": "2026-10-26T16:55:01.256057Z", + "start_time": "2026-10-26T16:55:01.237276Z" } }, "cell_type": "code", @@ -433,18 +565,16 @@ "id": "0f8fb861", "outputs": [ { - "ename": "NameError", - "evalue": "name 'abfe_default_settings' is not defined", - "output_type": "error", - "traceback": [ - "\u001B[31m---------------------------------------------------------------------------\u001B[39m", - "\u001B[31mNameError\u001B[39m Traceback (most recent call last)", - "\u001B[36mCell\u001B[39m\u001B[36m \u001B[39m\u001B[32mIn[1]\u001B[39m\u001B[32m, line 17\u001B[39m\n\u001B[32m 14\u001B[39m \u001B[38;5;28mprint\u001B[39m(\u001B[33mf\u001B[39m\u001B[33m\"\u001B[39m\u001B[33mTotal simulation time: \u001B[39m\u001B[38;5;132;01m{\u001B[39;00mabfe_total_sim_time.m\u001B[38;5;132;01m:\u001B[39;00m\u001B[33m.2f\u001B[39m\u001B[38;5;132;01m}\u001B[39;00m\u001B[33m ns\u001B[39m\u001B[33m\"\u001B[39m)\n\u001B[32m 15\u001B[39m \u001B[38;5;28;01mreturn\u001B[39;00m abfe_total_sim_time\n\u001B[32m---> \u001B[39m\u001B[32m17\u001B[39m default_abfe_cost = estimate_abfe_septop_cost(\u001B[43mabfe_default_settings\u001B[49m)\n", - "\u001B[31mNameError\u001B[39m: name 'abfe_default_settings' is not defined" + "name": "stdout", + "output_type": "stream", + "text": [ + "complex: 5.75 ns MD + 330.0 ns replica exchange = 335.75 ns\n", + "solvent: 0.8 ns MD + 154.0 ns replica exchange = 154.8 ns\n", + "Total simulation time: 1471.65 ns\n" ] } ], - "execution_count": 1 + "execution_count": 9 }, { "metadata": {}, @@ -459,8 +589,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.965069Z", - "start_time": "2026-10-26T11:45:16.913325Z" + "end_time": "2026-10-26T16:55:01.289159Z", + "start_time": "2026-10-26T16:55:01.268956Z" } }, "cell_type": "code", @@ -501,7 +631,7 @@ ] } ], - "execution_count": 9 + "execution_count": 10 }, { "metadata": {}, @@ -518,8 +648,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.218775Z", - "start_time": "2026-10-26T11:45:17.183233Z" + "end_time": "2026-10-26T16:55:01.322404Z", + "start_time": "2026-10-26T16:55:01.304106Z" } }, "cell_type": "code", @@ -559,7 +689,7 @@ ] } ], - "execution_count": 12 + "execution_count": 11 }, { "metadata": {}, @@ -574,8 +704,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.251783Z", - "start_time": "2026-10-26T11:45:17.234232Z" + "end_time": "2026-10-26T16:55:01.357944Z", + "start_time": "2026-10-26T16:55:01.325925Z" } }, "cell_type": "code", @@ -588,7 +718,114 @@ ], "id": "6c9de9ae4d246d9c", "outputs": [], - "execution_count": 13 + "execution_count": 12 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "#### Increasing the number of independent repeats\n", + "\n", + "The following example shows how to increase the number of independent repeats per single execution of the `AbsoluteBindingProtocol` to `4`. The number of repeats is shared by both legs, and every repeat runs the pre-equilibration MD simulation for each leg again, so the cost increases linearly with the number of repeats. As with the `RelativeHybridTopologyProtocol`, you can also execute the protocol multiple times with a single repeat per execution, see the [quickrun guide](https://docs.openfree.energy/en/stable/guide/execution/quickrun_execution.html#parallel-execution-of-repeats-with-quickrun) for more information." + ], + "id": "d1e2f3a4b5c6a01" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-27T09:21:28.528751Z", + "start_time": "2026-10-27T09:21:28.078100Z" + } + }, + "cell_type": "code", + "source": [ + "abfe_repeats_settings = AbsoluteBindingProtocol.default_settings()\n", + "\n", + "# Increase the number of independent repeats to 4\n", + "abfe_repeats_settings.protocol_repeats = 4\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "abfe_repeats_protocol = AbsoluteBindingProtocol(settings=abfe_repeats_settings)\n", + "\n", + "print(f\"Updated number of independent repeats: {abfe_repeats_protocol.settings.protocol_repeats}\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_abfe_cost = estimate_abfe_septop_cost(abfe_repeats_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_abfe_cost - default_abfe_cost) / default_abfe_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "d1e2f3a4b5c6a02", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated number of independent repeats: 4\n", + "complex: 5.75 ns MD + 330.0 ns replica exchange = 335.75 ns\n", + "solvent: 0.8 ns MD + 154.0 ns replica exchange = 154.8 ns\n", + "Total simulation time: 1962.20 ns\n", + "Percentage increase in total simulation time: 33.33%\n" + ] + } + ], + "execution_count": 20 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "#### Increasing the equilibration length\n", + "\n", + "The `AbsoluteBindingProtocol` has two equilibration stages per leg: the initial MD simulation controlled by `*_equil_simulation_settings`, which is used to select the restraint atoms, and the equilibration performed at each lambda window before production sampling, controlled by `*_simulation_settings.equilibration_length`. The following example doubles the equilibration length at each lambda window for both legs. As the length is set per lambda window, the cost increases linearly with the number of windows." + ], + "id": "d1e2f3a4b5c6a03" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-27T09:22:18.217432Z", + "start_time": "2026-10-27T09:22:18.145871Z" + } + }, + "cell_type": "code", + "source": [ + "abfe_equilibration_settings = AbsoluteBindingProtocol.default_settings()\n", + "\n", + "# Increase the equilibration length at each lambda window to 2 ns for both legs\n", + "abfe_equilibration_settings.complex_simulation_settings.equilibration_length = 2 * unit.nanoseconds\n", + "abfe_equilibration_settings.solvent_simulation_settings.equilibration_length = 2 * unit.nanoseconds\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "abfe_equilibration_protocol = AbsoluteBindingProtocol(settings=abfe_equilibration_settings)\n", + "\n", + "print(f\"Updated complex equilibration length: {abfe_equilibration_protocol.settings.complex_simulation_settings.equilibration_length.m} ns\")\n", + "print(f\"Updated solvent equilibration length: {abfe_equilibration_protocol.settings.solvent_simulation_settings.equilibration_length.m} ns\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_abfe_cost = estimate_abfe_septop_cost(abfe_equilibration_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_abfe_cost - default_abfe_cost) / default_abfe_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "d1e2f3a4b5c6a04", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated complex equilibration length: 2 ns\n", + "Updated solvent equilibration length: 2 ns\n", + "complex: 5.75 ns MD + 360.0 ns replica exchange = 365.75 ns\n", + "solvent: 0.8 ns MD + 168.0 ns replica exchange = 168.8 ns\n", + "Total simulation time: 1603.65 ns\n", + "Percentage increase in total simulation time: 8.97%\n" + ] + } + ], + "execution_count": 21 }, { "metadata": {}, @@ -603,8 +840,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.632171Z", - "start_time": "2026-10-26T11:45:16.525149Z" + "end_time": "2026-10-26T16:55:01.409452Z", + "start_time": "2026-10-26T16:55:01.359409Z" } }, "cell_type": "code", @@ -644,7 +881,7 @@ ] } ], - "execution_count": 4 + "execution_count": 13 }, { "metadata": {}, @@ -659,8 +896,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:16.850622Z", - "start_time": "2026-10-26T11:45:16.790688Z" + "end_time": "2026-10-27T09:22:33.285791Z", + "start_time": "2026-10-27T09:22:33.245289Z" } }, "cell_type": "code", @@ -677,7 +914,7 @@ ] } ], - "execution_count": 7 + "execution_count": 22 }, { "metadata": {}, @@ -692,8 +929,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.016519Z", - "start_time": "2026-10-26T11:45:16.983332Z" + "end_time": "2026-10-26T16:55:01.460079Z", + "start_time": "2026-10-26T16:55:01.438377Z" } }, "cell_type": "code", @@ -733,7 +970,7 @@ ] } ], - "execution_count": 10 + "execution_count": 15 }, { "metadata": {}, @@ -748,8 +985,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.286028Z", - "start_time": "2026-10-26T11:45:17.264318Z" + "end_time": "2026-10-27T09:22:52.631446Z", + "start_time": "2026-10-27T09:22:52.537454Z" } }, "cell_type": "code", @@ -790,13 +1027,13 @@ ] } ], - "execution_count": 14 + "execution_count": 23 }, { "metadata": { "ExecuteTime": { - "end_time": "2026-10-26T11:45:17.315038Z", - "start_time": "2026-10-26T11:45:17.307686Z" + "end_time": "2026-10-27T09:22:53.372160Z", + "start_time": "2026-10-27T09:22:53.350053Z" } }, "cell_type": "code", @@ -809,7 +1046,114 @@ ], "id": "32e3f17b3a1e849", "outputs": [], - "execution_count": 15 + "execution_count": 24 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "#### Increasing the number of independent repeats\n", + "\n", + "The following example shows how to increase the number of independent repeats per single execution of the `SepTopProtocol` to `4`. The number of repeats is shared by both legs, and every repeat runs the pre-equilibration MD simulation for each leg again, so the cost increases linearly with the number of repeats. As with the `RelativeHybridTopologyProtocol`, you can also execute the protocol multiple times with a single repeat per execution, see the [quickrun guide](https://docs.openfree.energy/en/stable/guide/execution/quickrun_execution.html#parallel-execution-of-repeats-with-quickrun) for more information." + ], + "id": "d1e2f3a4b5c6b01" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-27T09:23:15.668621Z", + "start_time": "2026-10-27T09:23:15.596538Z" + } + }, + "cell_type": "code", + "source": [ + "septop_repeats_settings = SepTopProtocol.default_settings()\n", + "\n", + "# Increase the number of independent repeats to 4\n", + "septop_repeats_settings.protocol_repeats = 4\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "septop_repeats_protocol = SepTopProtocol(settings=septop_repeats_settings)\n", + "\n", + "print(f\"Updated number of independent repeats: {septop_repeats_protocol.settings.protocol_repeats}\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_septop_cost = estimate_abfe_septop_cost(septop_repeats_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_septop_cost - default_septop_cost) / default_septop_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "d1e2f3a4b5c6b02", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated number of independent repeats: 4\n", + "complex: 2.2 ns MD + 209.0 ns replica exchange = 211.2 ns\n", + "solvent: 2.2 ns MD + 297.0 ns replica exchange = 299.2 ns\n", + "Total simulation time: 2041.60 ns\n", + "Percentage increase in total simulation time: 33.33%\n" + ] + } + ], + "execution_count": 25 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "#### Increasing the equilibration length\n", + "\n", + "The `SepTopProtocol` has two equilibration stages per leg: the initial MD simulation controlled by `*_equil_simulation_settings`, which is used to select the restraint atoms, and the equilibration performed at each lambda window before production sampling, controlled by `*_simulation_settings.equilibration_length`. The following example doubles the equilibration length at each lambda window for both legs. As the length is set per lambda window, the cost increases linearly with the number of windows." + ], + "id": "d1e2f3a4b5c6b03" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-10-27T09:23:20.742926Z", + "start_time": "2026-10-27T09:23:20.671284Z" + } + }, + "cell_type": "code", + "source": [ + "septop_equilibration_settings = SepTopProtocol.default_settings()\n", + "\n", + "# Increase the equilibration length at each lambda window to 2 ns for both legs\n", + "septop_equilibration_settings.complex_simulation_settings.equilibration_length = 2 * unit.nanoseconds\n", + "septop_equilibration_settings.solvent_simulation_settings.equilibration_length = 2 * unit.nanoseconds\n", + "\n", + "# Create a new protocol object with the updated settings\n", + "septop_equilibration_protocol = SepTopProtocol(settings=septop_equilibration_settings)\n", + "\n", + "print(f\"Updated complex equilibration length: {septop_equilibration_protocol.settings.complex_simulation_settings.equilibration_length.m} ns\")\n", + "print(f\"Updated solvent equilibration length: {septop_equilibration_protocol.settings.solvent_simulation_settings.equilibration_length.m} ns\")\n", + "\n", + "# Calculate the new total simulation time\n", + "new_septop_cost = estimate_abfe_septop_cost(septop_equilibration_settings)\n", + "\n", + "# Print the percentage increase in total simulation time\n", + "percentage_increase = ((new_septop_cost - default_septop_cost) / default_septop_cost) * 100\n", + "print(f\"Percentage increase in total simulation time: {percentage_increase.m:.2f}%\")" + ], + "id": "d1e2f3a4b5c6b04", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated complex equilibration length: 2 ns\n", + "Updated solvent equilibration length: 2 ns\n", + "complex: 2.2 ns MD + 228.0 ns replica exchange = 230.2 ns\n", + "solvent: 2.2 ns MD + 324.0 ns replica exchange = 326.2 ns\n", + "Total simulation time: 1669.20 ns\n", + "Percentage increase in total simulation time: 9.01%\n" + ] + } + ], + "execution_count": 26 } ], "metadata": {