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, and MonteCarloBarostat 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 as lambda 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 standard NonbondedForce 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 with lambda (and are covered by the NonbondedForce 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 separate CustomNonbondedForce that handles electrostatic interactions that change with lambda. This force uses interaction groups to exclude interactions that do not change with lambda (and are covered by the NonbondedForce 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 of perses.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 to True. 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 only lambda will be exposed. Otherwise, each individual energy term will have a lambda 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