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
auto_constraints_factor	The factor by which to multiply the timestep when deciding whether or not an auto-constrained bond should be constrained. The calculated period of a bond vibration is compared to the simulation timestep multiplied by this factor. If the period is less than this, the bond is constrained. This defaults to 10.
barostat_frequency	The frequency at which the barostat acts to perform the MC moves to change the box volume when performing constant pressure simulations (default 25).
check_for_h_by_ambertype	Boolean value, e.g. `True` or `False` as to whether hydrogen atoms can be detected based on a H? amber type. This is the default `True`.
check_for_h_by_element	Boolean value, e.g. `True` or `False` as to whether hydrogen atoms can be detected based on the element. This is the default `True`.
check_for_h_by_mass	Boolean value, e.g. `True` or `False` as to whether hydrogen atoms can be detected based on their mass. This is the default `True`.
check_for_h_by_max_mass	Boolean value, e.g. `True` or `False` as to whether hydrogen atoms can be detected based on the maximum mass of the end states. This is the default `True`.
com_reset_frequency	The frequency at which the `CMMotionRemover` acts to remove center of mass relative motion. If this is not set (the default) then center of mass motion is not removed.
constraint	Type of constraint to use for bonds and/or angles. Valid strings are `none`, `auto-bonds`, `h-bonds`, `h-bonds-not-perturbed`, `h-bonds-not-heavy-perturbed`, `bonds`, `bonds-not-perturbed`, `bonds-not-heavy-perturbed`, `h-bonds-h-angles`, `h-bonds-h-angles-not-perturbed`, `h-bonds-h-angles-not-heavy-perturbed`, `bonds-h-angles`, `bonds-h-angles-not-perturbed` and ``bonds-h-angles-not-heavy-perturbed
coulomb_power	The coulomb power parameter used by the softening potential used to soften electrostatic interactions involving ghost atoms. This defaults to 0.
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`
dynamic_constraints	Whether or not the constraints applied to perturbable bonds should be updated with λ (defaults to `True`).
fixed	The atoms in the system that should be fixed (not moved)
ignore_perturbations	Whether or not to ignore any perturbations and only set up a perturbable molecule as a non-perurbable molecule from only the reference state.
include_constrained_energies	Whether or not to include the bond and angle energies of bonds and angles that are constrained. This defaults to `True`, as we normally do want to calculate all of the energies of the internals, even if they are constrained.
integrator	The MD integrator to use, e.g. `verlet`, `leapfrog`, `langevin`, `langevin_middle`, `nose_hoover`, `brownian`, `andersen`
friction	Friction value for the integrator, in inverse time, e.g. `5.0 / sr.units.picosecond`
lambda	The λ-value at which to set up the system (assuming this contains any perturbable molecules or restraints)
perturbable_constraint	The constraint to use for perturbable molecules. These are the same options as `constraint`, and will override that choice for perturbable molecules if this is set.
platform	Any valid OpenMM platform string, e.g. `CUDA`, `OpenCL`, `Metal`, `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`
restraints	The `Restraints` object (or list of objects) that describe the restraints that should be added to the system.
schedule	The `LambdaSchedule` to use that controls how parameters are modified with λ
shift_coulomb	The coulomb delta parameter used by the softening potential used to soften electrostatic interactions involving ghost atoms. This defaults to 1.0 Å.
shift_delta	The shift_delta parameter to use for the softening potential used to soften LJ interactions involving ghost atoms. This defaults to 2.0 Å.
space	Space in which the simulation should be conducted, e.g. sr.vol.Cartesian
swap_end_states	Whether to swap the end states of a perturbable molecule (i.e. treat the perturbed state as the reference state and vice versa). This defaults to False.
taylor_power	The taylor power parameter used by the taylor algorithm for the softening potential used to soften LJ interactions involving ghost atoms. This defaults to 1.
temperature	Any temperature value, e.g. `25*sr.units.celsius`
threads	The number of threads to use in the CPU platform
timestep	Time between integration steps, e.g. `2 * sr.units.femtosecond`
tolerance	The tolerance to use for the PME calculation, e.g. `0.0001`
use_dispersion_correction	Whether or not to use the dispersion correction to deal with cutoff issues. This is very expensive.
use_taylor_softening	Whether or not to use the taylor algorithm to soften interactions involving ghost atoms. This defaults to False.
use_zacharias_softening	Whether or not to use the zacharias algorithm to soften interactions involving ghost atoms. This defaults to True. Note that one of zacharias or taylor softening must be True, with zacharias taking precedence.

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 most of the OpenMM options via arguments to the dynamics() function, e.g.

>>> d = mols.dynamics(temperature="25oC")
>>> print(d.ensemble())
canonical (NVT) ensemble { temperature = 298.15 C }

… note:

The function will automatically convert strings to units if these are
needed, e.g. ``25oC`` will automatically be converted to 25 Celsius.

>>> d = mols.dynamics(timestep="4fs", lambda_value=0.5)

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

>>> d = mols.dynamics(map={"temperature": "25oC"})
>>> 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(temperature="25oC", 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(timestep="1fs")
>>> print(d.constraint())
none
>>> d = mols.dynamics(timestep="5fs")
>>> 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 }
)

Running dynamics and saving frames and energies#

You can run dynamics via the run() function, e.g.

>>> d = mols.dynamics(timestep="4fs", temperature="25oC")
>>> d.run("100ps")

would run 100 picoseconds of dynamics.

At the end, you can extract the final system using the commit() function, e.g.

>>> mols = d.commit()

You can set the frequency at which trajectory frames and energies are saved via the save_frequency argument, e.g.

>>> d.run("100ps", save_frequency="10ps")

would save energies and trajectory frames every 10 picoseconds. You can specifiy different frequencies for these via the energy_frequency and/or frame_frequency arguments, e.g.

>>> d.run("1ns", energy_frequency="1ps", frame_frequency="100ps")

would save energies every picosecond and frames every 100 picoseconds.

By default, only coordinates are saved. You can choose to save velocities as well by setting save_velocities=True, e.g.

>>> d.run("10ps", save_frequency="1ps", save_velocities=True)

By default, energies are saved only for the simulated λ-value of the system. You can request energies to be saved for other λ-values using the lambda_windows argument, e.g.

>>> d.run("100ps", energy_frequency="1ps", lambda_windows=[0.0, 0.5, 1.0])

would save the energies at λ-values 0.0, 0.5 and 1.0 for every picosecond of the trajectory. You can pass in as many or few λ-windows as you wish.

The coordinate/velocity frames are saved to the trajectory property of the molecules, and are accessible identically to trajectories loaded from files (e.g. via that property of the .trajectory() function).

The energies are saved to the energy_trajectory property of the returned molecules, and accessible via that property or the energy_trajectory() function.