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 |
|
|---|---|
Module |
|
Settings |
|
MD Engine |
[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, the default settings can be printed clearly as a dictionary:
[4]:
settings.dict()
[4]:
{'forcefield_settings': {'constraints': 'hbonds',
'rigid_water': True,
'remove_com': False,
'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'},
'thermo_settings': {'temperature': 298.15 <Unit('kelvin')>,
'pressure': 0.9869232667160129 <Unit('standard_atmosphere')>,
'ph': None,
'redox_potential': None},
'system_settings': {'nonbonded_method': 'PME',
'nonbonded_cutoff': 1.0 <Unit('nanometer')>},
'solvation_settings': {'solvent_model': 'tip3p',
'solvent_padding': 1.2 <Unit('nanometer')>},
'alchemical_settings': {'lambda_functions': 'default',
'lambda_windows': 11,
'unsampled_endstates': False,
'use_dispersion_correction': False,
'softcore_LJ_v2': True,
'softcore_alpha': 0.85,
'interpolate_old_and_new_14s': False,
'flatten_torsions': False},
'alchemical_sampler_settings': {'online_analysis_interval': 250,
'n_repeats': 3,
'sampler_method': 'repex',
'online_analysis_target_error': 0.0 <Unit('boltzmann_constant * kelvin')>,
'online_analysis_minimum_iterations': 500,
'flatness_criteria': 'logZ-flatness',
'gamma0': 1.0,
'n_replicas': 11},
'engine_settings': {'compute_platform': None},
'integrator_settings': {'timestep': 4 <Unit('femtosecond')>,
'collision_rate': 1.0 <Unit('1 / picosecond')>,
'n_steps': 250 <Unit('timestep')>,
'reassign_velocities': False,
'n_restart_attempts': 20,
'constraint_tolerance': 1e-06,
'barostat_frequency': 25 <Unit('timestep')>},
'simulation_settings': {'equilibration_length': 1.0 <Unit('nanosecond')>,
'production_length': 5.0 <Unit('nanosecond')>,
'forcefield_cache': 'db.json',
'minimization_steps': 5000,
'output_filename': 'simulation.nc',
'output_structure': 'hybrid_system.pdb',
'output_indices': 'not water',
'checkpoint_interval': 250 <Unit('timestep')>,
'checkpoint_storage': 'checkpoint.nc'}}
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(
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
remove_com=False, # False: Do not remove center of mass motion
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',
),
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
),
system_settings=equil_rfe_settings.SystemSettings(
nonbonded_method='PME', # Particle Mesh Ewald for long range electrostatics
nonbonded_cutoff=1.0 * unit.nm, # Cut off Lennard-Jones interactions beyond 1 nm
),
solvation_settings=equil_rfe_settings.SolvationSettings(
solvent_model='tip3p', # Solvent model to generate starting coords
solvent_padding=1.2 * unit.nm, # Total distance between periodic image starting coords
),
alchemical_settings=equil_rfe_settings.AlchemicalSettings(
lambda_functions='default', # Interpolation functions for force field parameters
lambda_windows=11, # Split the transformation over n lambda windows
unsampled_endstates=False, # False: Use only the explicit lambda windows λ∈[0,-1]
use_dispersion_correction=False, # False: Use LJ dispersion correction only at endpoints
softcore_LJ_v2=True, # True: Use LJ potential from Gapsys et al. (JCTC 2012)
softcore_alpha=0.85, # Set soft-core Lennard-Jones potential parameter α
interpolate_old_and_new_14s=False, # False: Keep all exceptions (1,4 or otherwise) at all λ
flatten_torsions=False, # False: Keep all torsions at all lambda windows
),
alchemical_sampler_settings=equil_rfe_settings.AlchemicalSamplerSettings(
n_repeats=3, # Run n independent repeats of each transformation
# H-REX Sampling settings
sampler_method='repex', # Sample lambda with Hamiltonian Replica Exchange
n_replicas=11, # Number of HREX replicas with sampler_method='repex'
# # SAMS sampling settings
# sampler_method='sams', # Sample lambda with Self Adjusted Mixture Sampling
# flatness_criteria='logZ-flatness', # Criteria for asymptotically optimal sampling
# gamma0=1.0, # Initial adaptation rate w/ sampler_method='SAMS'
# Compute & write out free energies every n MCMC attempts:
online_analysis_interval=250,
# Compute & write out free energies only after the first n MCMC attempts:
online_analysis_minimum_iterations=500,
# Don't stop sampling early, no matter how low the estimated error gets:
online_analysis_target_error=0.0 * unit.boltzmann_constant * unit.kelvin
# # Stop sampling when estimated error is small enough:
# online_analysis_target_error=0.2 * unit.boltzmann_constant * unit.kelvin,
),
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
collision_rate=1.0 / unit.picosecond, # Langevin integration collision rate γ
n_steps=250 * unit.timestep, # Attempt an MCMC move every n integration steps
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
),
simulation_settings=equil_rfe_settings.SimulationSettings(
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
output_filename='simulation.nc', # Filename to save trajectory
output_structure='hybrid_system.pdb', # Filename to save starting coordinates
checkpoint_storage='checkpoint.nc', # 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.timestep, # Save a checkpoint every n integration steps
),
)
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.