diff --git a/cookbook/exploring_rbfe_settings.ipynb b/cookbook/exploring_rbfe_settings.ipynb new file mode 100644 index 0000000..212ad01 --- /dev/null +++ b/cookbook/exploring_rbfe_settings.ipynb @@ -0,0 +1,759 @@ +{ + "cells": [ + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "# 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", + "This section covers the settings that affect sampling, computational cost and statistical precision.\n", + "\n", + "### Inspecting the default settings\n", + "\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 a `Protocol` can be inspected via the common `default_settings()` class method. The following example shows how to retrieve the default settings for the `RelativeHybridTopologyProtocol`:" + ] + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-02T15:16:06.567836Z", + "start_time": "2026-11-02T15:15:50.439285Z" + } + }, + "cell_type": "code", + "source": [ + "from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol\n", + "\n", + "settings = RelativeHybridTopologyProtocol.default_settings()" + ], + "id": "5690f07a5cd27245", + "outputs": [], + "execution_count": 1 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "
\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 setting for the `RelativeHybridTopologyProtocol` are:\n", + "\n", + "| Setting | Description | Setting | Type |\n", + "|---------------------------------|--------------------------------------------------------------------------------------------------------------------------|--------------------------------------------|----------|\n", + "| **Lambda windows** | The number of discrete alchemical states simulated during the transformation. | `lambda_settings.lambda_windows` | `int` |\n", + "| **Production length** | The length of the production simulation performed at each lambda window. | `simulation_settings.production_length` | `openff.units.unit.Quantity` |\n", + "| **Equilibration length** | The length of the equilibration simulation performed at each lambda window using the selected multisate sampling method. | `simulation_settings.equilibration_length` | `openff.units.unit.Quantity` |\n", + "| **Protocol repeats** | The number of independent executions of the protocol. | `protocol_repeats` | `int` |\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:" + ], + "id": "94c914fc75cd8fcb" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-02T15:20:32.141958Z", + "start_time": "2026-11-02T15:20:32.071883Z" + } + }, + "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": 3 + }, + { + "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", + "This formula applies to all multistate protocols including:\n", + "\n", + "- `RelativeHybridTopologyProtocol`\n", + "- `AbsoluteBindingProtocol`\n", + "- `SepTopProtocol`\n", + "- `AbsoluteSolvationProtocol`\n", + "\n", + "So for the default settings the total simulation time is:" + ], + "id": "704d5dc6e6f895f5" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-02T16:31:19.687229Z", + "start_time": "2026-11-02T16:31:19.639276Z" + } + }, + "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", + "\n", + " total_sim_time = n_repeats * n_windows * (production_length + equilibration_length)\n", + " return total_sim_time\n", + "\n", + "total_sim_time_per_leg = estimate_rbfe_cost(settings)\n", + "print(f\"Total simulation time: {total_sim_time_per_leg}\")" + ], + "id": "fd81e13e2c34684d", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total simulation time: 198.0 nanosecond\n" + ] + } + ], + "execution_count": 18 + }, + { + "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" + }, + { + "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 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": "771ba6b9" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-02T15:21:27.259600Z", + "start_time": "2026-11-02T15:21:27.184581Z" + } + }, + "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_per_leg = estimate_rbfe_cost(protocol.settings)\n", + "print(f\"New total simulation time: {new_total_sim_time_per_leg}\")\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": "2cf9b3c7030df36f", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated production length: 10 nanosecond\n", + "New total simulation time: 363.0 nanosecond\n", + "Percentage increase in total simulation time: 83.33%\n" + ] + } + ], + "execution_count": 11 + }, + { + "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", + "
\n", + "Warning: Running additional lambda windows for a complete simulation is not supported. Increasing the number of lambda windows 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 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" + ], + "id": "3d6bb208a91ce6d0" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-02T15:21:38.243515Z", + "start_time": "2026-11-02T15:21:38.194029Z" + } + }, + "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", + "\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}\")\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": "a8a0a2b5ba644eb4", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated number of lambda windows: 22\n", + "New total simulation time: 396.0 nanosecond\n", + "Percentage increase in total simulation time: 100.00%\n" + ] + } + ], + "execution_count": 12 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Increasing the number of independent repeats\n", + "\n", + "\n", + "Independent repeats run the same protocol from the same input structures but different initial velocities. By default, **openfe** reports the **average** free energy change of the independent repeats as its aggregated estimate and the sample standard deviation is used to quantify the uncertainty of the prediction. Adding more repeats is a simple way to refine the estimate of the mean and uncertainty, but at a linear increase in computational cost. See the best practices paper by [May et al.](https://livecomsjournal.org/index.php/livecoms/article/view/v2i1e18378) for more information.\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-11-03T09:30:33.223560Z", + "start_time": "2026-11-03T09:30:32.220024Z" + } + }, + "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}\")\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 nanosecond\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 are far from equilibrium and the system would benefit from additional relaxation before production sampling begins.\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-11-02T15:22:03.407650Z", + "start_time": "2026-11-02T15:22:03.359829Z" + } + }, + "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}\")\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}\")\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 nanosecond\n", + "New total simulation time: 231.0 nanosecond\n", + "Percentage increase in total simulation time: 16.67%\n" + ] + } + ], + "execution_count": 15 + }, + { + "cell_type": "markdown", + "id": "c1d2e3f4a5b60001", + "metadata": {}, + "source": [ + "### Other protocols\n", + "\n", + "The above concepts apply to most free energy protocols, but the settings are located differently. Here we will briefly cover the differences in setting layout for the following protocols:\n", + "\n", + "- `AbsoluteBindingProtocol`\n", + "- `SepTopProtocol`\n", + "\n", + "These protocols execute both the `complex` and `solvent` legs within a single transformation, the settings to control the production sampling are then prefixed with the leg name `complex_*` or `solvent_*`.\n", + "\n", + "In addition, each leg has a `*_equil_simulation_settings` object controlling an initial MD simulation used to determine the restraint atoms in the complex and to seed both legs production simulations.\n" + ] + }, + { + "cell_type": "code", + "id": "other-protocol-imports", + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-03T10:23:12.015819Z", + "start_time": "2026-11-03T10:23:11.851672Z" + } + }, + "source": [ + "from openfe.protocols.openmm_afe import AbsoluteBindingProtocol\n", + "from openfe.protocols.openmm_septop import SepTopProtocol\n" + ], + "outputs": [], + "execution_count": 28 + }, + { + "cell_type": "markdown", + "id": "other-protocol-initial-md", + "metadata": {}, + "source": [ + "#### Initial MD equilibration settings\n", + "\n", + "Both protocols expose the initial MD settings for each leg as `complex_equil_simulation_settings` and `solvent_equil_simulation_settings`. These control the initial MD simulation (including NVT equilibration), not equilibration at each lambda window. The same setting names and approach apply to both protocols; here we update the NVT equilibration for both legs of an `AbsoluteBindingProtocol`.\n" + ] + }, + { + "cell_type": "code", + "id": "other-protocol-initial-md-example", + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-03T10:23:12.897324Z", + "start_time": "2026-11-03T10:23:12.835092Z" + } + }, + "source": [ + "abfe_settings = AbsoluteBindingProtocol.default_settings()\n", + "\n", + "for leg in (\"complex\", \"solvent\"):\n", + " equil_settings = getattr(abfe_settings, f\"{leg}_equil_simulation_settings\")\n", + " equil_settings.equilibration_length_nvt = 2 * unit.nanoseconds\n", + "\n", + "abfe_protocol = AbsoluteBindingProtocol(settings=abfe_settings)\n", + "print(f\"Updated complex initial NVT equilibration: {abfe_protocol.settings.complex_equil_simulation_settings.equilibration_length_nvt}\")\n", + "print(f\"Updated solvent initial NVT equilibration: {abfe_protocol.settings.solvent_equil_simulation_settings.equilibration_length_nvt}\")\n" + ], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Updated complex initial NVT equilibration: 2 nanosecond\n", + "Updated solvent initial NVT equilibration: 2 nanosecond\n" + ] + } + ], + "execution_count": 29 + }, + { + "cell_type": "markdown", + "id": "c1d2e3f4a5b60003", + "metadata": {}, + "source": [ + "#### Estimating computational cost\n", + "\n", + "For both protocols, each leg runs an initial MD simulation once per repeat, followed by equilibration and production at each lambda window. The same estimate applies to both:\n", + "\n", + "$$T_{total} = N_{repeats} \\times \\sum_{legs} \\left[ T_{MD} + N_{\\lambda windows} \\times (T_{production} + T_{equilibration}) \\right]$$\n", + "\n", + "Here, $T_{MD}$ is the total initial MD time set by `*_equil_simulation_settings`; the per-window lengths and window counts come from `*_simulation_settings` and the lambda schedule.\n" + ] + }, + { + "cell_type": "code", + "id": "other-protocol-cost-example", + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-03T10:23:13.847637Z", + "start_time": "2026-11-03T10:23:13.791142Z" + } + }, + "source": [ + "def estimate_abfe_septop_cost(settings):\n", + " 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", + " total_sim_time += md_time + replica_exchange_time\n", + "\n", + " return total_sim_time * settings.protocol_repeats\n", + "\n", + "abfe_default_settings = AbsoluteBindingProtocol.default_settings()\n", + "septop_default_settings = SepTopProtocol.default_settings()\n", + "print(f\"ABFE default total simulation time: {estimate_abfe_septop_cost(abfe_default_settings).m:.2f} ns\")\n", + "print(f\"SepTop default total simulation time: {estimate_abfe_septop_cost(septop_default_settings).m:.2f} ns\")\n" + ], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ABFE default total simulation time: 1471.65 ns\n", + "SepTop default total simulation time: 1531.20 ns\n" + ] + } + ], + "execution_count": 30 + }, + { + "cell_type": "markdown", + "id": "abfe-lambda-explanation", + "metadata": {}, + "source": [ + "### AbsoluteBindingProtocol - Increasing the number of lambda windows\n", + "\n", + "\n", + "ABFE uses one list of lambda values per scaled interaction (`lambda_elec`, `lambda_vdw` and `lambda_restraints`) allowing for exact control over the lambda spacing as well as number of windows. This requires you to define one value per window, with all lists for a leg having the same length, and that length must match `n_replicas` for the leg.\n" + ] + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-03T10:25:16.442991Z", + "start_time": "2026-11-03T10:25:16.417102Z" + } + }, + "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)" + ], + "id": "b1c2d3e4f5a60002", + "outputs": [], + "execution_count": 35 + }, + { + "cell_type": "markdown", + "id": "abfe-lambda-validation-tip", + "metadata": {}, + "source": [ + "
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 validate method to check that the lambda schedules are valid for all legs of the calculation, use them to validate your changes and catch issues before execution. This does require a representative system however, see the example below for more details.\n", + "
" + ] + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-03T10:25:17.577688Z", + "start_time": "2026-11-03T10:25:17.081423Z" + } + }, + "cell_type": "code", + "source": [ + "# Validate the lambda schedules for the protocol, this requires a valid inputs for the transformation\n", + "from rdkit import Chem\n", + "import openfe\n", + "\n", + "# Define a solvent for the calculation\n", + "solvent = openfe.SolventComponent()\n", + "# Load the protein and ligand from a pdb\n", + "protein = openfe.ProteinComponent.from_pdb_file(\"assets/tyk2_protein.pdb\")\n", + "# Load the ligands from an SDF file\n", + "supplier = Chem.SDMolSupplier(\"assets/tyk2_ligands.sdf\", removeHs=False)\n", + "ligands = [openfe.SmallMoleculeComponent.from_rdkit(mol) for mol in supplier if mol is not None]\n", + "\n", + "# Build an example chemical system to validate the lambda schedule\n", + "system_a = openfe.ChemicalSystem({\n", + " \"ligand\": ligands[0],\n", + " \"protein\": protein,\n", + " \"solvent\": solvent\n", + "})\n", + "\n", + "system_b = openfe.ChemicalSystem({\n", + " \"protein\": protein,\n", + " \"solvent\": solvent\n", + "})\n", + "\n", + "abfe_protocol = AbsoluteBindingProtocol(settings=abfe_window_settings)\n", + "abfe_protocol.validate(\n", + " stateA=system_a,\n", + " stateB=system_b,\n", + " mapping=None, # The mapping is not required for ABFE or SepTop, but is required for RBFE\n", + ")" + ], + "id": "6c9de9ae4d246d9c", + "outputs": [], + "execution_count": 36 + }, + { + "cell_type": "markdown", + "id": "septop-lambda-explanation", + "metadata": {}, + "source": [ + "### SepTopProtocol - Increasing the number of lambda windows\n", + "\n", + "SepTop defines separate lambda lists for each scaled interaction and for each of the states A and B: `lambda_elec_A/B`, `lambda_vdw_A/B` and `lambda_restraints_A/B`. All lists for a leg must have the same length, and the matching `n_replicas` must equal the number of windows.\n" + ] + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-03T10:23:19.692203Z", + "start_time": "2026-11-03T10:23:19.665574Z" + } + }, + "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" + ], + "id": "b1c2d3e4f5a60004", + "outputs": [], + "execution_count": 33 + }, + { + "cell_type": "markdown", + "id": "septop-lambda-validation-tip", + "metadata": {}, + "source": "Validate the edited schedule before running the protocol using the `validate` method of the protocol, this requires a representative system to be provided.\n" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-11-03T10:47:35.479549Z", + "start_time": "2026-11-03T10:47:35.085448Z" + } + }, + "cell_type": "code", + "source": [ + "# Validate the updated lambda schedule for the protocol using an example system with the inputs from the ABFE cell above\n", + "system_a = openfe.ChemicalSystem({\n", + " \"ligand\": ligands[0],\n", + " \"protein\": protein,\n", + " \"solvent\": solvent\n", + "})\n", + "system_b = openfe.ChemicalSystem({\n", + " \"ligand\": ligands[1],\n", + " \"protein\": protein,\n", + " \"solvent\": solvent\n", + "})\n", + "\n", + "septop_protocol = SepTopProtocol(settings=septop_window_settings)\n", + "septop_protocol.validate(\n", + " stateA=system_a,\n", + " stateB=system_b,\n", + " mapping=None, # The mapping is not required for ABFE or SepTop, but is required for RBFE\n", + ")" + ], + "id": "32e3f17b3a1e849", + "outputs": [], + "execution_count": 38 + } + ], + "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 +}