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.
from openff.units import unit
Protocol
¶Your choice of Protocol
determines how free energy sampling is performed. Here, we will be looking into the RelativeHybridTopologyProtocol
:
Name | RelativeHybridTopologyProtocol |
---|---|
Module | openfe.protocols.openmm_rfe |
Settings | RelativeHybridTopologyProtocolSettings |
MD Engine | OpenMM |
from openfe.protocols.openmm_rfe import (
RelativeHybridTopologyProtocol,
RelativeHybridTopologyProtocolSettings,
)
from openfe.protocols.openmm_rfe import equil_rfe_settings
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:
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:
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.1.1'}, '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:
settings.simulation_settings.production_length = 10.0 * unit.nanosecond
Alternatively, settings can be specified by hand when creating the settings object:
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.1.1',
# 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
),
)
# Double check that the above settings match the defaults - delete this cell if you configure things yourself!
# See https://github.com/OpenFreeEnergy/ExampleNotebooks/issues/138
#assert settings == RelativeHybridTopologyProtocol.default_settings()
--------------------------------------------------------------------------- AssertionError Traceback (most recent call last) Cell In[15], line 3 1 # Double check that the above settings match the defaults - delete this cell if you configure things yourself! 2 # See https://github.com/OpenFreeEnergy/ExampleNotebooks/issues/138 ----> 3 assert settings == RelativeHybridTopologyProtocol.default_settings() AssertionError:
default_settings = RelativeHybridTopologyProtocol.default_settings()
for name, s in settings:
print(name)
if getattr(settings, name) != getattr(default_settings, name):
print(f"{settings}, {name} != \n {default_settings}, {name}")
forcefield_settings 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=<Quantity(1.0, 'nanometer')>, nonbonded_method='PME') thermo_settings=ThermoSettings(temperature=<Quantity(298.15, 'kelvin')>, pressure=<Quantity(0.986923267, 'standard_atmosphere')>, ph=None, redox_potential=None) protocol_repeats=3 solvation_settings=OpenMMSolvationSettings(solvent_model='tip3p', solvent_padding=<Quantity(1.2, 'nanometer')>) 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=<Quantity(0.8, 'nanometer')>, 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=<Quantity(1.0, 'nanosecond')>, production_length=<Quantity(5.0, 'nanosecond')>, minimization_steps=5000, time_per_iteration=<Quantity(1.0, 'picosecond')>, real_time_analysis_interval=<Quantity(250.0, 'picosecond')>, early_termination_target_error=<Quantity(0.0, 'kilocalorie / mole')>, real_time_analysis_minimum_time=<Quantity(500.0, 'picosecond')>, 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=<Quantity(4.0, 'femtosecond')>, langevin_collision_rate=<Quantity(1.0, '1 / picosecond')>, barostat_frequency=<Quantity(25.0, 'timestep')>, remove_com=False, reassign_velocities=False, n_restart_attempts=20, constraint_tolerance=1e-06) output_settings=MultiStateOutputSettings(checkpoint_interval=<Quantity(250.0, 'picosecond')>, forcefield_cache='db.json', output_indices='not water', checkpoint_storage_filename='checkpoint.chk', output_filename='simulation.nc', output_structure='hybrid_system.pdb'), forcefield_settings != 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=<Quantity(1.0, 'nanometer')>, nonbonded_method='PME') thermo_settings=ThermoSettings(temperature=<Quantity(298.15, 'kelvin')>, pressure=<Quantity(0.986923267, 'standard_atmosphere')>, ph=None, redox_potential=None) protocol_repeats=3 solvation_settings=OpenMMSolvationSettings(solvent_model='tip3p', solvent_padding=<Quantity(1.2, 'nanometer')>) 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=<Quantity(0.8, 'nanometer')>, 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=<Quantity(1.0, 'nanosecond')>, production_length=<Quantity(5.0, 'nanosecond')>, minimization_steps=5000, time_per_iteration=<Quantity(1, 'picosecond')>, real_time_analysis_interval=<Quantity(250, 'picosecond')>, early_termination_target_error=<Quantity(0.0, 'kilocalorie_per_mole')>, real_time_analysis_minimum_time=<Quantity(500, 'picosecond')>, 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=<Quantity(4, 'femtosecond')>, langevin_collision_rate=<Quantity(1.0, '1 / picosecond')>, barostat_frequency=<Quantity(25, 'timestep')>, remove_com=False, reassign_velocities=False, n_restart_attempts=20, constraint_tolerance=1e-06) output_settings=MultiStateOutputSettings(checkpoint_interval=<Quantity(250, 'picosecond')>, forcefield_cache='db.json', output_indices='not water', checkpoint_storage_filename='checkpoint.chk', output_filename='simulation.nc', output_structure='hybrid_system.pdb'), forcefield_settings thermo_settings protocol_repeats solvation_settings partial_charge_settings lambda_settings alchemical_settings simulation_settings engine_settings integrator_settings output_settings
Protocol
¶However you produce the ProtocolSettings
object, the final Protocol
can be constructed from the settings object:
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.