Alchemical transformations¶
Tools for nonequilibrium alchemical transformations
Nonequilibrium switching¶
Relative alchemical transformations¶
Principle of HybridTopologyFactory¶
HybridTopologyFactory is a class that automates the construction of so-called hybrid topologies and systems. In short, this amounts to using an atom map to merge
two initial systems into a new system that contains a union of the former systems’ degrees of freedom, along with a lambda
parameter to control the degree to which the hybrid represents the old or new system. The process of creating this system happens in several steps, detailed below. Before proceeding, there are several important caveats:
- Atom maps:
Atoms whose constraints change between new and old system may not be mapped
Virtual sites are only permitted if their parameters are identical in the old and new systems
- Systems:
Only
HarmonicBondForce
,HarmonicAngleForce
,PeriodicTorsionForce
,NonbondedForce
, andMonteCarloBarostat
are supported. The presence of other forces will raise an exception.
Basic Terminology¶
Throughout the rest of this document, the following terms will be referenced:
Hybrid topology: An object that “fuses” two systems based on an atom map, and whose identity is controlled by a set of global parameters called
lambda
. Mapped atoms are switched from the type in one system to the type in the other; unmapped atoms are switched on or off, depending on which system they belong to.Old system: We use this to refer to the endpoint represented by
lambda=0
.New system: We use this to refer to the endpoint represented by
lambda=1
.
Assignment of Particles to Appropriate Groups¶
In order to properly assign parameters to the particles in the hybrid system, the particles are first assigned to one of four groups:
- Environment: These atoms are mapped between the new and old systems, and additionally are not part of a residue that differs between new and old systems.
Examples: solvent molecules, protein residues outside a changing residue
- Core: These atoms are mapped between the new and old system, but are part of a residue that differs between new and old systems. As such, they may need to change parameters.
Example: The carbon atoms in a benzene that is being transformed to a chlorobenzene.
- Unique old: These atoms are not in the map, and are present only in the old system
Example: an extraneous hydrogen atom in a cyclohexane being transformed into a benzene
- Unique new: These atoms are not in the map, and are present only in the new system
Example: The chlorine atom that appears when a benzene is transformed into a chlorobenzene.
Treatment of different interaction groups¶
For each supported force, we have to create at least one custom force, as well as an unmodified force. In general, the interactions are handled as follows:
Environment-Environment: These interactions are handled by regular unmodified forces, as their parameters never change.
Environment-{Core, Unique old, Unique new}: These interactions are handled by Custom*Forces. For interactions involving the core atoms, the interaction scales from the parameter in the old system to the parameter in the new system as
lambda
is switched from 0 to 1. For interactions involving the unique atoms, the interaction scales from full to zero aslambda
switches from 0 to 1 in unique old atoms, and vice versa for unique new atoms.Core-Core interactions: These interactions scale from old system parameters to new system parameters as
lambda
is switched from 0 to 1.Core-{unique old, unique new} interactions: These interactions scale from the parameters in the old system to the parameters in the new system as
lambda
is switched from 0 to 1. This means that unique atoms always maintain intramolecular interactions.Unique-old and unique-new interactions: These are disabled.
intra-Unique atoms: These interactions are always enabled, and as such are handled by an unmodified force.
Force terms¶
Given the above division of atoms into different groups, the assignment of interactions to different force classes is straightforward.
Bonds¶
Bonds are handled with a HarmonicBondForce
and a CustomBondForce
.
HarmonicBondForce
: This handles all bond interactions that do not change with lambda.CustomBondForce
: This handles bond interactions that change with lambda. Bond terms never become zero; they are always on to ensure that molecules do not fly apart
Angles¶
Angles are also handled with a HarmonicAngleForce
and a CustomAngleForce
HarmonicAngleForce
: This handles all angle interactions that do not change with lambda.CustomAngleForce
: This handles angle interactions that change with lambda. Angle terms never become zero; they are always on to ensure overlap.
Torsions¶
Torsions are handled with a PeriodicTorsionForce
and a CustomTorsionForce
PeriodicTorsionForce
: This handles all torsion interactions that do not change with lambda.CustomTorsionForce
: This handles interactions that change with lambda. Although the torsion potential for a given group of atoms is never set to zero, the force assignments sometimes include zero as an endpoint. This is to handle the case of torsion multiplicities greater than one.
Nonbonded (both sterics and electrostatics)¶
NonbondedForce
: This handles sterics and electrostatics interactions that do not change. All interactions that are involved in custom nonbonded forces are excluded from this force.CustomBondForce
: This handles the exceptions from the standardNonbondedForce
above.
Nonbonded (Sterics)¶
CustomNonbondedForce
: This handles softcore sterics for interactions that change with lambda. This force uses interaction groups to exclude interactions that do not change withlambda
(and are covered by theNonbondedForce
above). This force supports no cutoff, reaction field, and PME. An important note relates to the dispersion correction; by default, it is not enabled to improve performance. However, it can be enabled as a constructor argument in the HybridTopologyFactory.
Nonbonded (Electrostatics)¶
CustomNonbondedForce
: This is a separateCustomNonbondedForce
that handles electrostatic interactions that change with lambda. This force uses interaction groups to exclude interactions that do not change withlambda
(and are covered by theNonbondedForce
above). This force supports the options nocutoff, reaction field, and PME.
Specific unsupported forces¶
Any force that is not in one of the above-mentioned classes will cause an exception to be raised. Some important examples of this are:
Amoeba forces
GB forces
CMMotionRemover
Additional Features¶
In addition to creating a hybrid topology dependent on a set of lambda_force
functions, the HybridTopologyFactory can also take advantage of OpenMM’s automatic differentiation to provide dU/dlambda. In order to use this
functionality, you must provide alchemical functions to the constructor of HybridTopologyFactory. These functions are described below.
Basic Usage of HybridTopologyFactory¶
Basic Construction of Factory¶
In order to use the HybridTopologyFactory, several objects are necessary:
perses.rjmc.topology_proposal.TopologyProposal
is a container class which holds so-called new and old systems, topologies, and respective atom maps. It can be created manually using the constructor, or can be obtained through a subcalss ofperses.rjmc.topology_proposal.ProposalEngine
.You need to supply current positions and new positions for the old and new systems, respectively. It is assumed that the atoms which are mapped will have the same positions. If you only have positions for the old system, you can generate (decent) positions for the new system using
perses.rjmc.geometry.FFAllAngleGeometryEngine
. You should minimize the hybrid topology that results after this to avoid clashes, however.The dispersion correct is calculated for the sterics force if
use_dispersion_correction
is set toTrue
. It is false by default; whether to enable this depends on what sort of calculation you are doing:Your calculation requires frequent changes to
lambda
parameters, such as nonequilibrium switching: Set the dispersion correction to false.Your calculation runs a simulation at a set of constant
lambda
values (that don’t change during the simulation), such as thermodynamic integration: Set the dispersion correction to True.
The alchemical functions argument can be supplied to the
HybridTopologyFactory
, in which case onlylambda
will be exposed. Otherwise, each individual energy term will have alambda
associated with it.
Alchemical functions¶
The functions argument allows you to specify the mathematical relationship between a main lambda value, and each individual lambda. Adding the functions here also automatically enables the computation of the derivative of the energy with respect to lambda
, useful for various tasks.
Some important notes about specifying the alchemical functions:
All
lambda
values must be specified. More specifically, the dictionary must have the following form:
functions_dictionary = {
'lambda_bonds' : 'f_bond(lambda)'
'lambda_angles' : 'f_angle(lambda)'
'lambda_torsions' : 'f_torsion(lambda)'
'lambda_sterics' : 'f_sterics(lambda)'
'lambda_electrostatics' : 'f_electrostatics(lambda)'
}
where f_{energy_name}()
is an arbitrary function of lambda
compliant with OpenMM Lepton syntax. See the OpenMM documentation for more information on what kinds of functions can be implemented there.
Example script¶
Below is a very simple example script that constructs a HybridTopologyFactory
object that can transform a benzene into catechol in vacuum.
#import needed functionality
from topology_proposal import SmallMoleculeSetProposalEngine, TopologyProposal
from perses.annihilation.relative import HybridTopologyFactory
import simtk.openmm.app as app
from openmoltools import forcefield_generators
#these are utility functions to rapidly create test systems.
from perses.utils.openeye import createOEMolFromIUPAC, createSystemFromIUPAC
from perses.utils.data import get_data_filename
#We'll generate the systems by IUPAC name
mol_name = "benzene"
ref_mol_name = "catechol"
#create the benzene molecule and system, as well as a molecule of catechol
m, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(mol_name)
refmol = iupac_to_oemol(ref_mol_name)
#The SmallMoleculeSetProposalEngine class needs smiles strings
initial_smiles = oechem.OEMolToSmiles(m)
final_smiles = oechem.OEMolToSmiles(refmol)
#set up the SystemGenerator
gaff_xml_filename = get_data_filename("data/gaff.xml")
system_generator = SystemGenerator([gaff_filename])
#The GeometryEngine will be useful for generating new positions
geometry_engine = FFAllAngleGeometryEngine()
#Instantiate the proposal engine, which will generate the TopologyProposal.
proposal_engine = SmallMoleculeSetProposalEngine(
[initial_smiles, final_smiles], system_generator, residue_name=mol_name)
#generate a TopologyProposal with the SmallMoleculeSetProposalEngine
topology_proposal = proposal_engine.propose(solvated_system, solvated_topology)
#Generate new positions with the GeometryEngine--note that the second return value (logp) is ignored for this purpose
new_positions, _ = geometry_engine.propose(topology_proposal, solvated_positions, beta)
#instantiate the HybridTopologyFactory
factory = HybridTopologyFactory(topology_proposal, old_positions, new_positions)
#Now you have the objects you need to run a simulation, controlled by the `lambda` parameters.
hybrid_system = factory.hybrid_system
hybrid_topology = factory.hybrid_topology
initial_hybrid_positions = factory.hybrid_positions