Examples using SMIRNOFF with the toolkit

The following examples are available in the openforcefield toolkit repository:

Index of provided examples

  • SMIRNOFF_simulation/ - simulation of a molecule in the gas phase with the SMIRNOFF forcefield format
  • SMIRNOFF_comparison/ - compare molecule energies from SMIRNOFF/OpenMM and AMBER
  • forcefield_modification/ - Jupyter notebook example of modifying a forcefield parameter and evaluating how it changes an energy
  • chemicalEnvironments/ - example and documentation of using chemical environment objects to manipulate environment being considered, generate example SMIRKS, etc. Also contains IPython notebook using the chemical environment for depiction.
  • mixedFF_structure - example of how to use SMIRNOFF parameters for small molecules in combination with more conventional force fields for proteins and other components of your system (using ParmEd to combine parameterized structures)
  • partial_bondorder/ - example of using partial bond orders for parameterizing molecules, and showing how the SMIRNOFF format extends; see https://github.com/openforcefield/smarty/issues/53
  • host_guest_simulation/ - example of setting up a basic MD simulation for a host-guest system using the SMIRNOFF force fields.
  • label_molecule/ - example showing how to check which parameters are used in which molecules, do parameter usage statistics, etc.

These examples are described in more detail below.

Simulate a molecule in the gas phase using SMIRNOFF forcefield format

This example code (run_molecule.py) demonstrates how to simulate a molecule in the gas phase using the SMIRNOFF forcefield format.

The essential idea is to create an openforcefield Topology object (as well as positions, if present) from a mol2 file:

# Load molecule
from openforcefield.topology import Molecule
molecule = Molecule.from_file(mol_filename)

# Create an openforcefield Topology
from openforcefield.topology import Topology
topology = Topology.from_molecules(molecule)

# Get positions in OpenMM-compatible format
positions = molecule.positions

We then load a SMIRNOFF forcefield and create an OpenMM System object to simulate:

# Load a SMIRNOFF small molecule forcefield for alkanes, ethers, and alcohols
forcefield = ForceField(offxml_filename)

# Create the OpenMM system
system = forcefield.create_system(topology)

Compare small molecule parm@frosst energies between SMIRNOFF and AMBER prmtop/inpcrd

This example demonstrates how to compare parm@frosst energies of molecules from the AlkEthOH set between systems generated from a SMIRNOFF .offxml file and systems generated from AMBER prmtop/inpcrd files.

The key idea here is to use the utility compare_molecule_energies to compare the molecule parameterized by ForceField with a corresponding system constructed via an AMBER prmtop/inpcrd pair:

# Load molecule
from openmmtools.topology import Molecule
molecule = Molecule.from_file(mol_filename)

# Load forcefield
from openforcefield.typing.engines.smirnoff import ForceField
forcefield = ForceField(offxml_filename)

# Compare energies
from openforcefield.tests.utils import compare_molecule_energies
results = compare_molecule_energies(prmtop_filename, inpcrd_filename, forcefield, molecule)

Manipulating SMIRNOFF parameters after loading into a ForceField object

In this example, we illustrate how to use the API for manipulating SMIRNOFF parameters after they are loaded into a ForceField object.

Combining a SMIRNOFF parameterized small molecule with an AMBER parameterized protein using ParmEd

This example illustrates how the ParmEd utility can be used to merge a small molecule parameterized by SMIRNOFF with a traditionally parameterized protein (or other biopolymer) to create a fully parameterized protein-ligand system.

The essential idea is to first create a ParmEd structure from the SMIRNOFF parameterized small molecule

# Load the small molecule
from openforcefield.utils import get_data_filename
ligand_filename = get_data_filename('molecules/toluene.mol2')
molecule = Molecule.from_file(ligand_filename)

# Load the smirnoff99Frosst force field
from openforcefield.typing.engines import smirnoff
forcefield = smirnoff.ForceField('smirnoff99Frosst.offxml')

# Create a ParmEd structure for the molecule
molecule_structure = forcefield.create_parmed_structure(topology=molecule.to_topology(), positions=molecule.positions)

Next, a ParmEd structure is created for the protein (parameterized with AMBER using the OpenMM app.ForceField implementation):

# Load the protein topology
from openforcefield.utils import generate_protein_structure
protein_pdb_filename = get_data_filename('proteins/T4-protein.pdb')
protein_pdbfile = app.PDBFile(protein_pdb_filename)

# Load the AMBER protein force field, along with a solvent force field
from simtk.openmm import app
protein_forcefield = 'amber99sbildn.xml'
solvent_forcefield = 'tip3p.xml'
forcefield = app.ForceField(protein_forcefield, solvent_forcefield)

# Parameterize the protein
# TODO: Can we create an OpenMM ffxml file for the ligand?
protein_system = forcefield.createSystem(proteinpdb.topology)

# Create a ParmEd Structure for the protein
protein_structure = parmed.openmm.load_topology(proteinpdb.topology,
                                                protein_system,
                                                xyz=proteinpdb.positions)

Finally, the two ParmEd Structure objects can be combined by concatenation:

# Combine the ParmEd Structure objects to produce a fully parameterized complex
complex_structure = protein_structure + molecule_structure

Using partial bond orders

In this example, we demonstrate how a SMIRNOFF force field can use partial bond order information to assign bond parameters relevant for benzene as well as singly and doubly-bonded trivalent carbon, all with a single bond entry.

The file Frosst_AlkEthOH_extracarbons.offxml is a SMIRNOFF FFXML file adding extra carbon parameters to cover benzene (adapted from work by Christopher Bayly) and uses partial bond orders to compress the [#6X3]-[#6X3] , [#6X3]:[#6X3] , and [#6X3]=[#6X3] bond parameters to a single line.

The <Bonds> section contains

<Bonds potential="harmonic" length_unit="angstroms" k_unit="kilocalories_per_mole/angstrom**2" fractional_bondorder_method="Wiberg" fractional_bondorder_interpolation="linear">
   ...
   <Bond smirks="[#6X3:1]!#[#6X3:2]" k_bondorder1="820.0" k_bondorder2="1098" length_bondorder1="1.45" length_bondorder2="1.35"/> <!-- Christopher Bayly from parm99, Aug 2016 -->
   <Bond smirks="[#6X3:1]-[#1:2]" k="734.0" length="1.080"/> <!-- Christopher Bayly from parm99, Aug 2016 -->
</Bonds>

The jupyter notebook test_partialbondorder.ipynb illustrates the application of these parameters to benzene, and inspection of the resulting bonded parameters.

Setting up a host-guest simulation using SMIRNOFF parameters

This example illustrates how easy it is to set up rather nontrivial simulations with SMIRNOFF force fields. In this case, we will begin with a guest SMILES string and a 3D structure of a host, dock the host into the guest, solvate, and proceed to simulations.

Authorship

This example was provided by David Mobley (UCI)

Source materials

Example files

Inspecting which SMIRNOFF parameters are assigned to specific molecules

This example illustrates how to inspect which SMIRNOFF parameters are assigned to specific atoms for a test molecule, using the AlkEthOH parameter set.

Example files

  • label_molecule.py: Creates a simple example molecule and uses the labeler to indicate what force terms are applied to which atoms
  • get_parameter_statistics.py: Pulls details on which parameters are used in each molecule, and which molecules each parameter occurs in