Choose and Configure a Protocol#

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.

Setup#

[1]:
from openff.units import unit

Choose a Protocol#

Your choice of Protocol determines how free energy sampling is performed. At present, the only protocol available is RelativeHybridTopologyProtocol:

Name

RelativeHybridTopologyProtocol

Module

openfe.protocols.openmm_rfe

Settings

RelativeHybridTopologyProtocolSettings

MD Engine

OpenMM

[2]:
from openfe.protocols.openmm_rfe import (
    RelativeHybridTopologyProtocol,
    RelativeHybridTopologyProtocolSettings,
)
from openfe.protocols.openmm_rfe import equil_rfe_settings

Configure Protocol Settings#

From the Defaults#

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:

[3]:
settings = RelativeHybridTopologyProtocol.default_settings()

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:

[4]:
settings
{'alchemical_settings': {'endstate_dispersion_correction': False,
                         'explicit_charge_correction': False,
                         'explicit_charge_correction_cutoff': <Quantity(0.8, 'nanometer')>,
                         'softcore_LJ': 'gapsys',
                         'softcore_alpha': 0.85,
                         'turn_off_core_unique_exceptions': False,
                         'use_dispersion_correction': False},
 'engine_settings': {'compute_platform': None},
 'forcefield_settings': {'constraints': 'hbonds',
                         'forcefields': ['amber/ff14SB.xml',
                                         'amber/tip3p_standard.xml',
                                         'amber/tip3p_HFE_multivalent.xml',
                                         'amber/phosaa10.xml'],
                         'hydrogen_mass': 3.0,
                         'nonbonded_cutoff': <Quantity(1.0, 'nanometer')>,
                         'nonbonded_method': 'PME',
                         'rigid_water': True,
                         'small_molecule_forcefield': 'openff-2.0.0'},
 'integrator_settings': {'barostat_frequency': <Quantity(25, 'timestep')>,
                         'constraint_tolerance': 1e-06,
                         'langevin_collision_rate': <Quantity(1.0, '1 / picosecond')>,
                         'n_restart_attempts': 20,
                         'reassign_velocities': False,
                         'remove_com': False,
                         'timestep': <Quantity(4, 'femtosecond')>},
 'lambda_settings': {'lambda_functions': 'default', 'lambda_windows': 11},
 'output_settings': {'checkpoint_interval': <Quantity(250, 'picosecond')>,
                     'checkpoint_storage_filename': 'checkpoint.chk',
                     'forcefield_cache': 'db.json',
                     'output_filename': 'simulation.nc',
                     'output_indices': 'not water',
                     'output_structure': 'hybrid_system.pdb'},
 'partial_charge_settings': {'nagl_model': None,
                             'number_of_conformers': None,
                             'off_toolkit_backend': 'ambertools',
                             'partial_charge_method': 'am1bcc'},
 'protocol_repeats': 3,
 'simulation_settings': {'early_termination_target_error': <Quantity(0.0, 'kilocalorie_per_mole')>,
                         'equilibration_length': <Quantity(1.0, 'nanosecond')>,
                         'minimization_steps': 5000,
                         'n_replicas': 11,
                         'production_length': <Quantity(5.0, 'nanosecond')>,
                         'real_time_analysis_interval': <Quantity(250, 'picosecond')>,
                         'real_time_analysis_minimum_time': <Quantity(500, 'picosecond')>,
                         'sampler_method': 'repex',
                         'sams_flatness_criteria': 'logZ-flatness',
                         'sams_gamma0': 1.0,
                         'time_per_iteration': <Quantity(1, 'picosecond')>},
 'solvation_settings': {'solvent_model': 'tip3p',
                        'solvent_padding': <Quantity(1.2, 'nanometer')>},
 'thermo_settings': {'ph': None,
                     'pressure': <Quantity(0.986923267, 'standard_atmosphere')>,
                     'redox_potential': None,
                     'temperature': <Quantity(298.15, 'kelvin')>}}

The production simulations could be lengthened from 5 ns to 10 ns:

[5]:
settings.simulation_settings.production_length = 10.0 * unit.nanosecond

From Scratch#

Alternatively, settings can be specified by hand when creating the settings object:

[6]:
settings = RelativeHybridTopologyProtocolSettings(
    protocol_repeats=3,  # Number of independent repeats of the Protocol transformation
    forcefield_settings=equil_rfe_settings.OpenMMSystemGeneratorFFSettings(
        constraints='hbonds',              # 'hbonds': Use constraints for bonds involving hydrogen
        rigid_water=True,                  # True: Use constraints for bonds in water
        hydrogen_mass=3.0,                 # Perform hydrogen mass repartitioning
        forcefields=[                      # OpenMM force fields to use for solvents and proteins
            'amber/ff14SB.xml',
            'amber/tip3p_standard.xml',
            'amber/tip3p_HFE_multivalent.xml',
            'amber/phosaa10.xml'
        ],

        # Small molecule force field to use with OpenMM template generator:
        small_molecule_forcefield='openff-2.0.0',

        # Nonbonded settings
        nonbonded_method='PME',            # Particle Mesh Ewald for long range electrostatics
        nonbonded_cutoff=1.0 * unit.nm,    # Cut off Lennard-Jones interactions beyond 1 nm
    ),
    thermo_settings=equil_rfe_settings.ThermoSettings(
        temperature=298.15 * unit.kelvin,  # Set thermostat temperature
        pressure=1 * unit.bar,             # Set barostat pressure
        ph=None,                           # None: Do not keep pH constant
        redox_potential=None               # None: Do not keep redox potential constant
    ),
    solvation_settings=equil_rfe_settings.OpenMMSolvationSettings(
        solvent_model='tip3p',             # Solvent model to generate starting coords
        solvent_padding=1.2 * unit.nm,     # Total distance between periodic image starting coords
    ),
    partial_charge_settings=equil_rfe_settings.OpenFFPartialChargeSettings(
        partial_charge_method='am1bcc',    # Partial charge method applied - am1bcc
        off_toolkit_backend='ambertools',  # Toolkit to use for partial charge assignment - ambertools
        number_of_conformers=None,         # None: use input conformer for partial charge assignment
        nagl_model=None,                   # None: not using NAGL so no model needs to be chosen
    ),
    lambda_settings=equil_rfe_settings.LambdaSettings(
        lambda_functions='default',        # Interpolation functions for force field parameters
        lambda_windows=11,                 # Split the transformation over n lambda windows
    ),
    alchemical_settings=equil_rfe_settings.AlchemicalSettings(
        # False: Don't use unsampled (non-hybrid) endstates for long range correction
        endstate_dispersion_correction=False,
        use_dispersion_correction=False,   # False: Don't use dispersion correction
        softcore_LJ='gapsys',              # Use LJ potential from Gapsys et al. (JCTC 2012)
        softcore_alpha=0.85,               # Set soft-core Lennard-Jones potential parameter α
        # False: Keep all exceptions (1,4 or otherwise) at all λ
        tun_off_core_unique_exceptions=False,

        # Explicit charge correction settings
        # False: don't apply explicit charge correction using an alchemical water
        explicit_charge_correction=False,
        # Cutoff distance for choosing alchemical waters
        explicit_charge_correction_cutoff=0.8 * unit.nm,
    ),
    simulation_settings=equil_rfe_settings.MultiStateSimulationSettings(
        # Simulation lengths
        minimization_steps=5000,                    # Minimize potential energy for n steps
        equilibration_length=1.0 * unit.nanosecond, # Simulation time to equilibrate for
        production_length=5.0 * unit.nanosecond,    # Simulation time to collect data for

        # Alchemical Space Sampling settings
        n_replicas=11,                           # Number of replicas sampling alchemical space
        sampler_method='repex',                  # Sample lambda with Hamiltonian Replica Exchange
        time_per_iteration=1*unit.ps,            # Time interval between state sampling (MCMC) attempts

        # SAMS sampling settings (used if sampler_method='sams')
        sams_flatness_criteria='logZ-flatness',  # Criteria for switch to asymptomatically optimal scheme
        sams_gamma0=1.0,                          # Initial SAMS weight adoption rate.

        # Settings to control free energy analysis
        # Time interval at which to perform an analysis of the free energies
        real_time_analysis_interval=250*unit.picosecond,
        # Minimum simulation time before energy analysis is carried out
        real_time_analysis_minimum_time=500*unit.picosecond,
        # Stop simulation if this target error is reached:
        early_termination_target_error=0.0*unit.kilocalorie_per_mole,
    ),
    engine_settings=equil_rfe_settings.OpenMMEngineSettings(
        compute_platform=None,              # Let OpenMM choose the best platform for your hardware
    ),
    integrator_settings=equil_rfe_settings.IntegratorSettings(
        timestep=4 * unit.femtosecond,         # Integration timestep
        langevin_collision_rate=1.0 / unit.picosecond,  # Langevin integrator collision rate γ
        reassign_velocities=False,             # False: Velocities are not lost through MCMC moves
        n_restart_attempts=20,                 # Restart simulations the first n times they blow up
        constraint_tolerance=1e-06,            # Tolerance for holonomic constraints
        barostat_frequency=25 * unit.timestep, # Attempt MC volume scaling every n integration steps
        remove_com=False,                      # False: don't remove the center of mass motion
    ),
    output_settings=equil_rfe_settings.MultiStateOutputSettings(
        output_filename='simulation.nc',            # Filename to save trajectory
        output_structure='hybrid_system.pdb',       # Filename to save starting coordinates
        checkpoint_storage_filename='checkpoint.chk',  # Filename for simulation checkpoints
        forcefield_cache='db.json',                 # Cache for small molecule residue templates
        output_indices='not water',                 # Do not save water positions
        checkpoint_interval=250 * unit.ps,          # Save a checkpoint every 250 picoseconds
    ),
)

Construct the Protocol#

However you produce the ProtocolSettings object, the final Protocol can be constructed from the settings object:

[8]:
protocol = RelativeHybridTopologyProtocol(settings)

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.