{ "cells": [ { "cell_type": "markdown", "id": "8efbded0-0b69-4392-b8df-7afbd817d733", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# Choose and Configure a Protocol" ] }, { "cell_type": "markdown", "id": "66156fbb-30e8-4d19-b1b0-1e488042ad22", "metadata": {}, "source": [ "A `Protocol` describes the simulation and sampling strategy for a free energy campaign. It is specified with subclasses of the `Protocol` class, and their associated `ProtocolSettings` subclasses." ] }, { "cell_type": "markdown", "id": "bbe95b7c-5396-4e0d-b1d3-21910ef2fb07", "metadata": {}, "source": [ "## Setup" ] }, { "cell_type": "code", "execution_count": 9, "id": "38b41f0a-d2a3-4231-bf08-1504d3824978", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "from openff.units import unit" ] }, { "cell_type": "markdown", "id": "893b1b18-3682-435f-a08f-e08c0f2f8cef", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Choose a `Protocol`" ] }, { "cell_type": "markdown", "id": "dd66b0da-7366-44d5-b02a-ba2a67581b55", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Your choice of `Protocol` determines how free energy sampling is performed. Here, we will be looking into the `RelativeHybridTopologyProtocol`:\n", "\n", "| Name | [RelativeHybridTopologyProtocol] |\n", "|:-------------|:-----------------------------------------|\n", "| Module | `openfe.protocols.openmm_rfe` |\n", "| Settings | [RelativeHybridTopologyProtocolSettings] | \n", "| MD Engine | [OpenMM] |\n", "\n", "[RelativeHybridTopologyProtocol]: https://docs.openfree.energy/en/stable/reference/api/generated/openfe.protocols.openmm_rfe.RelativeHybridTopologyProtocol.html\n", "[RelativeHybridTopologyProtocolSettings]: https://docs.openfree.energy/en/stable/reference/api/openmm_rfe.html#openfe.protocols.openmm_rfe.RelativeHybridTopologyProtocolSettings\n", "[OpenMM]: https://openmm.org" ] }, { "cell_type": "code", "execution_count": 10, "id": "b0757f5f-e38f-43ac-b43f-ee2bbea95f0d", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "from openfe.protocols.openmm_rfe import (\n", " RelativeHybridTopologyProtocol,\n", " RelativeHybridTopologyProtocolSettings,\n", ")\n", "from openfe.protocols.openmm_rfe import equil_rfe_settings" ] }, { "cell_type": "markdown", "id": "f9502961-84ce-4983-9c4f-9e07eb619d1f", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Configure Protocol Settings" ] }, { "cell_type": "markdown", "id": "39517e1a-0e03-414e-8475-59bfb96f3203", "metadata": {}, "source": [ "### From the Defaults" ] }, { "cell_type": "markdown", "id": "03fa0695-7066-494a-9e9d-6f7d7f3953d5", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The user-configurable settings of a `Protocol` are stored in a separate object that inherits from `ProtocolSettings`. The default settings object for a protocol can be retrieved with the `Protocol.default_settings` class method:" ] }, { "cell_type": "code", "execution_count": 11, "id": "96c32ead-a8d0-4283-950b-5d27c74ebf16", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "settings = RelativeHybridTopologyProtocol.default_settings()" ] }, { "cell_type": "markdown", "id": "b85b39a6-effd-423e-b8a5-b86bc2e7ee77", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The settings object is a Pydantic data class, and so can be edited and inspected in the usual ways. For example, you can call the object to print its contents as a dictionary:" ] }, { "cell_type": "code", "execution_count": 12, "id": "8d5a6281-50f9-4062-8fa4-375ce2df1bd6", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'alchemical_settings': {'endstate_dispersion_correction': False,\n", " 'explicit_charge_correction': False,\n", " 'explicit_charge_correction_cutoff': ,\n", " 'softcore_LJ': 'gapsys',\n", " 'softcore_alpha': 0.85,\n", " 'turn_off_core_unique_exceptions': False,\n", " 'use_dispersion_correction': False},\n", " 'engine_settings': {'compute_platform': None},\n", " 'forcefield_settings': {'constraints': 'hbonds',\n", " 'forcefields': ['amber/ff14SB.xml',\n", " 'amber/tip3p_standard.xml',\n", " 'amber/tip3p_HFE_multivalent.xml',\n", " 'amber/phosaa10.xml'],\n", " 'hydrogen_mass': 3.0,\n", " 'nonbonded_cutoff': ,\n", " 'nonbonded_method': 'PME',\n", " 'rigid_water': True,\n", " 'small_molecule_forcefield': 'openff-2.1.1'},\n", " 'integrator_settings': {'barostat_frequency': ,\n", " 'constraint_tolerance': 1e-06,\n", " 'langevin_collision_rate': ,\n", " 'n_restart_attempts': 20,\n", " 'reassign_velocities': False,\n", " 'remove_com': False,\n", " 'timestep': },\n", " 'lambda_settings': {'lambda_functions': 'default', 'lambda_windows': 11},\n", " 'output_settings': {'checkpoint_interval': ,\n", " 'checkpoint_storage_filename': 'checkpoint.chk',\n", " 'forcefield_cache': 'db.json',\n", " 'output_filename': 'simulation.nc',\n", " 'output_indices': 'not water',\n", " 'output_structure': 'hybrid_system.pdb'},\n", " 'partial_charge_settings': {'nagl_model': None,\n", " 'number_of_conformers': None,\n", " 'off_toolkit_backend': 'ambertools',\n", " 'partial_charge_method': 'am1bcc'},\n", " 'protocol_repeats': 3,\n", " 'simulation_settings': {'early_termination_target_error': ,\n", " 'equilibration_length': ,\n", " 'minimization_steps': 5000,\n", " 'n_replicas': 11,\n", " 'production_length': ,\n", " 'real_time_analysis_interval': ,\n", " 'real_time_analysis_minimum_time': ,\n", " 'sampler_method': 'repex',\n", " 'sams_flatness_criteria': 'logZ-flatness',\n", " 'sams_gamma0': 1.0,\n", " 'time_per_iteration': },\n", " 'solvation_settings': {'solvent_model': 'tip3p',\n", " 'solvent_padding': },\n", " 'thermo_settings': {'ph': None,\n", " 'pressure': ,\n", " 'redox_potential': None,\n", " 'temperature': }}\n" ] } ], "source": [ "settings" ] }, { "cell_type": "markdown", "id": "b39aec57-d263-473e-a314-ce2e4519bef9", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The production simulations could be lengthened from 5 ns to 10 ns:" ] }, { "cell_type": "code", "execution_count": 13, "id": "2de6c6a9-bceb-4c9c-ae5a-3d88a72c801a", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "settings.simulation_settings.production_length = 10.0 * unit.nanosecond" ] }, { "cell_type": "markdown", "id": "49cf72a9-5230-4232-8a04-001aa725a9b8", "metadata": {}, "source": [ "### From Scratch" ] }, { "cell_type": "markdown", "id": "606b0be8-6d76-4a5f-b9b0-9b54967ca988", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Alternatively, settings can be specified by hand when creating the settings object:" ] }, { "cell_type": "code", "execution_count": 14, "id": "844e6fa5-f363-47fa-881d-00328faa601c", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "settings = RelativeHybridTopologyProtocolSettings(\n", " protocol_repeats=3, # Number of independent repeats of the Protocol transformation\n", " forcefield_settings=equil_rfe_settings.OpenMMSystemGeneratorFFSettings(\n", " constraints='hbonds', # 'hbonds': Use constraints for bonds involving hydrogen\n", " rigid_water=True, # True: Use constraints for bonds in water\n", " hydrogen_mass=3.0, # Perform hydrogen mass repartitioning\n", " forcefields=[ # OpenMM force fields to use for solvents and proteins\n", " 'amber/ff14SB.xml',\n", " 'amber/tip3p_standard.xml',\n", " 'amber/tip3p_HFE_multivalent.xml',\n", " 'amber/phosaa10.xml'\n", " ],\n", " \n", " # Small molecule force field to use with OpenMM template generator:\n", " small_molecule_forcefield='openff-2.1.1',\n", " \n", " # Nonbonded settings\n", " nonbonded_method='PME', # Particle Mesh Ewald for long range electrostatics\n", " nonbonded_cutoff=1.0 * unit.nm, # Cut off Lennard-Jones interactions beyond 1 nm\n", " ),\n", " thermo_settings=equil_rfe_settings.ThermoSettings(\n", " temperature=298.15 * unit.kelvin, # Set thermostat temperature\n", " pressure=1 * unit.bar, # Set barostat pressure\n", " ph=None, # None: Do not keep pH constant\n", " redox_potential=None # None: Do not keep redox potential constant\n", " ),\n", " solvation_settings=equil_rfe_settings.OpenMMSolvationSettings(\n", " solvent_model='tip3p', # Solvent model to generate starting coords\n", " solvent_padding=1.2 * unit.nm, # Total distance between periodic image starting coords\n", " ),\n", " partial_charge_settings=equil_rfe_settings.OpenFFPartialChargeSettings(\n", " partial_charge_method='am1bcc', # Partial charge method applied - am1bcc\n", " off_toolkit_backend='ambertools', # Toolkit to use for partial charge assignment - ambertools\n", " number_of_conformers=None, # None: use input conformer for partial charge assignment\n", " nagl_model=None, # None: not using NAGL so no model needs to be chosen\n", " ),\n", " lambda_settings=equil_rfe_settings.LambdaSettings(\n", " lambda_functions='default', # Interpolation functions for force field parameters\n", " lambda_windows=11, # Split the transformation over n lambda windows\n", " ),\n", " alchemical_settings=equil_rfe_settings.AlchemicalSettings(\n", " # False: Don't use unsampled (non-hybrid) endstates for long range correction\n", " endstate_dispersion_correction=False,\n", " use_dispersion_correction=False, # False: Don't use dispersion correction\n", " softcore_LJ='gapsys', # Use LJ potential from Gapsys et al. (JCTC 2012)\n", " softcore_alpha=0.85, # Set soft-core Lennard-Jones potential parameter α\n", " # False: Keep all exceptions (1,4 or otherwise) at all λ\n", " tun_off_core_unique_exceptions=False,\n", " \n", " # Explicit charge correction settings\n", " # False: don't apply explicit charge correction using an alchemical water\n", " explicit_charge_correction=False,\n", " # Cutoff distance for choosing alchemical waters\n", " explicit_charge_correction_cutoff=0.8 * unit.nm,\n", " ),\n", " simulation_settings=equil_rfe_settings.MultiStateSimulationSettings(\n", " # Simulation lengths\n", " minimization_steps=5000, # Minimize potential energy for n steps\n", " equilibration_length=1.0 * unit.nanosecond, # Simulation time to equilibrate for\n", " production_length=5.0 * unit.nanosecond, # Simulation time to collect data for\n", "\n", " # Alchemical Space Sampling settings\n", " n_replicas=11, # Number of replicas sampling alchemical space\n", " sampler_method='repex', # Sample lambda with Hamiltonian Replica Exchange\n", " time_per_iteration=1*unit.ps, # Time interval between state sampling (MCMC) attempts\n", " \n", " # SAMS sampling settings (used if sampler_method='sams')\n", " sams_flatness_criteria='logZ-flatness', # Criteria for switch to asymptomatically optimal scheme\n", " sams_gamma0=1.0, # Initial SAMS weight adoption rate.\n", " \n", " # Settings to control free energy analysis\n", " # Time interval at which to perform an analysis of the free energies\n", " real_time_analysis_interval=250*unit.picosecond,\n", " # Minimum simulation time before energy analysis is carried out\n", " real_time_analysis_minimum_time=500*unit.picosecond,\n", " # Stop simulation if this target error is reached:\n", " early_termination_target_error=0.0*unit.kilocalorie_per_mole,\n", " ),\n", " engine_settings=equil_rfe_settings.OpenMMEngineSettings(\n", " compute_platform=None, # Let OpenMM choose the best platform for your hardware\n", " ),\n", " integrator_settings=equil_rfe_settings.IntegratorSettings(\n", " timestep=4 * unit.femtosecond, # Integration timestep\n", " langevin_collision_rate=1.0 / unit.picosecond, # Langevin integrator collision rate γ\n", " reassign_velocities=False, # False: Velocities are not lost through MCMC moves\n", " n_restart_attempts=20, # Restart simulations the first n times they blow up\n", " constraint_tolerance=1e-06, # Tolerance for holonomic constraints\n", " barostat_frequency=25 * unit.timestep, # Attempt MC volume scaling every n integration steps\n", " remove_com=False, # False: don't remove the center of mass motion\n", " ),\n", " output_settings=equil_rfe_settings.MultiStateOutputSettings(\n", " output_filename='simulation.nc', # Filename to save trajectory\n", " output_structure='hybrid_system.pdb', # Filename to save starting coordinates\n", " checkpoint_storage_filename='checkpoint.chk', # Filename for simulation checkpoints\n", " forcefield_cache='db.json', # Cache for small molecule residue templates\n", " output_indices='not water', # Do not save water positions\n", " checkpoint_interval=250 * unit.ps, # Save a checkpoint every 250 picoseconds\n", " ),\n", ")" ] }, { "cell_type": "code", "execution_count": 15, "id": "b499d852-4d7e-47c5-9342-31017baed37c", "metadata": { "nbsphinx": "hidden", "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "ename": "AssertionError", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[15], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;66;03m# Double check that the above settings match the defaults - delete this cell if you configure things yourself!\u001b[39;00m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;66;03m# See https://github.com/OpenFreeEnergy/ExampleNotebooks/issues/138\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m settings \u001b[38;5;241m==\u001b[39m RelativeHybridTopologyProtocol\u001b[38;5;241m.\u001b[39mdefault_settings()\n", "\u001b[0;31mAssertionError\u001b[0m: " ] } ], "source": [ "# Double check that the above settings match the defaults - delete this cell if you configure things yourself!\n", "# See https://github.com/OpenFreeEnergy/ExampleNotebooks/issues/138\n", "#assert settings == RelativeHybridTopologyProtocol.default_settings()" ] }, { "cell_type": "code", "execution_count": 17, "id": "3a6d3f17-2049-4456-a315-ce40ebd8969c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "forcefield_settings\n", "forcefield_settings=OpenMMSystemGeneratorFFSettings(constraints='hbonds', rigid_water=True, hydrogen_mass=3.0, forcefields=['amber/ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml', 'amber/phosaa10.xml'], small_molecule_forcefield='openff-2.0.0', nonbonded_cutoff=, nonbonded_method='PME') thermo_settings=ThermoSettings(temperature=, pressure=, ph=None, redox_potential=None) protocol_repeats=3 solvation_settings=OpenMMSolvationSettings(solvent_model='tip3p', solvent_padding=) partial_charge_settings=OpenFFPartialChargeSettings(partial_charge_method='am1bcc', off_toolkit_backend='ambertools', number_of_conformers=None, nagl_model=None) lambda_settings=LambdaSettings(lambda_functions='default', lambda_windows=11) alchemical_settings=AlchemicalSettings(softcore_LJ='gapsys', explicit_charge_correction_cutoff=, endstate_dispersion_correction=False, use_dispersion_correction=False, softcore_alpha=0.85, turn_off_core_unique_exceptions=False, explicit_charge_correction=False) simulation_settings=MultiStateSimulationSettings(equilibration_length=, production_length=, minimization_steps=5000, time_per_iteration=, real_time_analysis_interval=, early_termination_target_error=, real_time_analysis_minimum_time=, sampler_method='repex', sams_flatness_criteria='logZ-flatness', sams_gamma0=1.0, n_replicas=11) engine_settings=OpenMMEngineSettings(compute_platform=None) integrator_settings=IntegratorSettings(timestep=, langevin_collision_rate=, barostat_frequency=, remove_com=False, reassign_velocities=False, n_restart_attempts=20, constraint_tolerance=1e-06) output_settings=MultiStateOutputSettings(checkpoint_interval=, forcefield_cache='db.json', output_indices='not water', checkpoint_storage_filename='checkpoint.chk', output_filename='simulation.nc', output_structure='hybrid_system.pdb'), forcefield_settings != \n", " forcefield_settings=OpenMMSystemGeneratorFFSettings(constraints='hbonds', rigid_water=True, hydrogen_mass=3.0, forcefields=['amber/ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml', 'amber/phosaa10.xml'], small_molecule_forcefield='openff-2.1.1', nonbonded_cutoff=, nonbonded_method='PME') thermo_settings=ThermoSettings(temperature=, pressure=, ph=None, redox_potential=None) protocol_repeats=3 solvation_settings=OpenMMSolvationSettings(solvent_model='tip3p', solvent_padding=) partial_charge_settings=OpenFFPartialChargeSettings(partial_charge_method='am1bcc', off_toolkit_backend='ambertools', number_of_conformers=None, nagl_model=None) lambda_settings=LambdaSettings(lambda_functions='default', lambda_windows=11) alchemical_settings=AlchemicalSettings(softcore_LJ='gapsys', explicit_charge_correction_cutoff=, endstate_dispersion_correction=False, use_dispersion_correction=False, softcore_alpha=0.85, turn_off_core_unique_exceptions=False, explicit_charge_correction=False) simulation_settings=MultiStateSimulationSettings(equilibration_length=, production_length=, minimization_steps=5000, time_per_iteration=, real_time_analysis_interval=, early_termination_target_error=, real_time_analysis_minimum_time=, sampler_method='repex', sams_flatness_criteria='logZ-flatness', sams_gamma0=1.0, n_replicas=11) engine_settings=OpenMMEngineSettings(compute_platform=None) integrator_settings=IntegratorSettings(timestep=, langevin_collision_rate=, barostat_frequency=, remove_com=False, reassign_velocities=False, n_restart_attempts=20, constraint_tolerance=1e-06) output_settings=MultiStateOutputSettings(checkpoint_interval=, forcefield_cache='db.json', output_indices='not water', checkpoint_storage_filename='checkpoint.chk', output_filename='simulation.nc', output_structure='hybrid_system.pdb'), forcefield_settings\n", "thermo_settings\n", "protocol_repeats\n", "solvation_settings\n", "partial_charge_settings\n", "lambda_settings\n", "alchemical_settings\n", "simulation_settings\n", "engine_settings\n", "integrator_settings\n", "output_settings\n" ] } ], "source": [ "default_settings = RelativeHybridTopologyProtocol.default_settings()\n", "\n", "for name, s in settings:\n", " print(name)\n", "\n", " if getattr(settings, name) != getattr(default_settings, name):\n", " print(f\"{settings}, {name} != \\n {default_settings}, {name}\")" ] }, { "cell_type": "markdown", "id": "99a37707-7658-4ccd-844e-a9a4cbe058aa", "metadata": {}, "source": [ "## Construct the `Protocol`" ] }, { "cell_type": "markdown", "id": "c9b77e53-67cc-4fc1-9386-e73cd005eb9b", "metadata": {}, "source": [ "However you produce the `ProtocolSettings` object, the final `Protocol` can be constructed from the settings object:" ] }, { "cell_type": "code", "execution_count": null, "id": "32b9b58a-d7ea-4071-9046-53e6990a1faa", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "protocol = RelativeHybridTopologyProtocol(settings)" ] }, { "cell_type": "markdown", "id": "cb02efa6-0949-4745-8bf7-b16580299ac3", "metadata": {}, "source": [ "Unlike `ProtocolSettings`, a `Protocol` instance is immutable. The only way to safely change the settings of a `Protocol` is to recreate it from the modified `ProtocolSettings` object." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.14" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }