OpenMM Integration#

The dynamics functionality in sire is based on tight integration with OpenMM.

This is achieved via code in sire.convert that can convert sire objects to their OpenMM equivalents.

You can convert any molecule, molecule view, collection of molecules or system objects into an OpenMM object using sire.convert.to() and passing in openmm as the desired format.

>>> import sire as sr
>>> mols = sr.load(sr.expand(sr.tutorial_url, "kigaki.gro", "kigaki.top"),
...                show_warnings=False)
>>> omm = sr.convert.to(mols, "openmm")
>>> print(omm)
<openmm.openmm.Context; proxy of <Swig Object of type 'OpenMM::Context *' at 0x2846e5cb0> >

The result is an OpenMM Context which can be used just like any other Context object.

For example, you can extract the OpenMM System via

>>> print(omm.getSystem())
<openmm.openmm.System; proxy of <Swig Object of type 'OpenMM::System *' at 0x2846e5d40> >

or the OpenMM Integrator

>>> print(omm.getIntegrator())
<openmm.openmm.VerletIntegrator; proxy of <Swig Object of type 'OpenMM::VerletIntegrator *' at 0x2846e5a70> >

or the OpenMM Platform

>>> print(omm.getPlatform())
<openmm.openmm.Platform; proxy of <Swig Object of type 'OpenMM::Platform *' at 0x2846e5ce0> >

as you could for any other OpenMM Context.

You could then recombine these with your own choice of Integrator or Platform if you want to use something different to what was created by sire.

Choosing options#

sire made a series of choices for various OpenMM options and parameters when it created the OpenMM Context. These choices were based on sensible defaults derived from what it could find within the molecules/system being converted.

You can override these choices by passing in a dictionary of key-value pairs as the map option to sire.convert.to().

For example, you can choose the integrator by setting a value for the integrator key, and set the temperature and pressure by setting values for the temperature and pressure keys.

>>> m = {"integrator": "langevin_middle",
...      "temperature": 25*sr.units.celsius,
...      "pressure": 1*sr.units.atm}
>>> omm = sr.convert.to(mols, "openmm", map=m)
>>> print(omm.getIntegrator())
<openmm.openmm.LangevinMiddleIntegrator; proxy of <Swig Object of type 'OpenMM::LangevinMiddleIntegrator *' at 0x295a07cc0> >

Available keys and allowable values are listed below.

Key

Valid values

constraint

Type of constraint to use for bonds and/or angles. Valid strings are none, h-bonds, bonds, h-bonds-h-angles and bonds-h-angles.

cutoff

Size of the non-bonded cutoff, e.g. 7.5*sr.units.angstrom

cutoff_type

Type of cutoff, e.g. PARTICLE_MESH_EWALD, PME, NO_CUTOFF, REACTION_FIELD, RF, EWALD

cpu_pme

Boolean value, e.g. True or False as to whether or not PME should be evaluated on the CPU rather than on the GPU.

device

Any valid OpenMM device number or device string, e.g. 0 would select device 0

dielectric

Dielectric value if a reaction field cutoff is used, e.g. 78.0

integrator

verlet, leapfrog, langevin, langevin_middle, nose_hoover, brownian

friction

Friction value for the integrator, in inverse time, e.g. 5.0 / sr.units.picosecond

platform

Any valid OpenMM platform string, e.g. CUDA, OpenCL, CPU, Reference

precision

Any valid OpenMM platform precision value, e.g. single, mixed or double.

pressure

Any pressure value, e.g. 1*sr.units.atm

space

Space in which the simulation should be conducted, e.g. sr.vol.Cartesian

temperature

Any temperature value, e.g. 25*sr.units.celsius

timestep

Time between integration steps, e.g. 2 * sr.units.femtosecond

Higher level API#

The Dynamics object and dynamics() function provides a higher level API for running molecular dynamics using the OpenMM Context created by sire.

You create a Dynamics object by calling the dynamics() function on the molecule, molecule view, collection or system that you want to simulate. For example

>>> d = mols.dynamics()

You can use this object to query the options that were passed into OpenMM.

>>> print(d.ensemble())
microcanonical (NVE) ensemble

You can set OpenMM options by passing the dictionary of key-value pairs as the map option.

>>> d = mols.dynamics(map={"temperature": 25*sr.units.celsius})
>>> print(d.ensemble())
canonical (NVT) ensemble { temperature = 298.15 C }

Note

sire automatically chooses the right OpenMM Integrator and barostat options based on the ensemble parameters.

It is a mistake to use an OpenMM Integrator that is not suited for the chosen ensemble.

>>> d = mols.dynamics(map={"temperature": 25*sr.units.celsius,
...                        "integrator": "verlet"})
ValueError: You cannot use a verlet integrator with the ensemble canonical (NVT) ensemble { temperature = 298.15 C }

You can also query other parameters.

>>> d = mols.dynamics(map={"timestep": 1*sr.units.femtosecond})
>>> print(d.constraint())
none
>>> d = mols.dynamics(map={"timestep": 5*sr.units.femtosecond})
>>> print(d.constraint())
bonds-h-angles
>>> print(d.timestep())
0.005 ps

Forcefield properties are automatically set based on the properties contained by the molecules. You can get a summary of these properties using the info() function.

>>> print(d.info())
ForceFieldInfo(
  space=PeriodicBox( ( 48.3263, 48.3263, 48.3263 ) ),
  cutoff_type=PME,
  cutoff=7.5 Å,
  params=Properties( tolerance => 0.0001 ),
  detail=MM ForceField{ amber::ff,
               combining_rules = arithmetic,
               1-4 scaling = 0.833333, 0.5,
               nonbonded = coulomb, lj,
               bond = harmonic, angle = harmonic,
               dihedral = cosine }
)

Some of these properties, such as detail, come from the forcefield parameters of the converted molecules. Others, such as the cutoff_type and cutoff are passed from the options given by the user (or derived as defaults). The space property is extracted from the System if that is passed, or is found from the space property from the first molecule that contains such a property. Sometimes, particularly if you aren’t using a System, it can be a good idea to manually set the space, e.g. to Cartesian if you are running a gas-phase simulation. In this case setting the cutoff_type to NO_CUTOFF will set the cutoff to a sufficiently large value so that the effect is that there is no cutoff. Setting the space to Cartesian will require disabling PME, as this cutoff type requires a periodic space. Instead, choose a cutoff type like reaction field.

>>> d = mols.dynamics(map={"space": sr.vol.Cartesian(),
...                        "cutoff_type": "NO_CUTOFF"})
>>> print(d.info())
ForceFieldInfo(
  space=Infinite cartesian space,
  cutoff_type=NO_CUTOFF,
  detail=MM ForceField{ amber::ff,
               combining_rules = arithmetic,
               1-4 scaling = 0.833333, 0.5,
               nonbonded = coulomb, lj,
               bond = harmonic, angle = harmonic,
               dihedral = cosine }
)
>>> d = mols.dynamics(map={"space": sr.vol.Cartesian(),
...                        "cutoff_type": "RF"})
>>> print(d.info())
ForceFieldInfo(
  space=Infinite cartesian space,
  cutoff_type=REACTION_FIELD,
  cutoff=7.5 Å,
  params=Properties( dielectric => 78.3 ),
  detail=MM ForceField{ amber::ff,
               combining_rules = arithmetic,
               1-4 scaling = 0.833333, 0.5,
               nonbonded = coulomb, lj,
               bond = harmonic, angle = harmonic,
               dihedral = cosine }
)