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 |
cutoff |
Size of the non-bonded cutoff, e.g.
|
cutoff_type |
Type of cutoff, e.g. |
cpu_pme |
Boolean value, e.g. |
device |
Any valid OpenMM device number or device string, e.g.
|
dielectric |
Dielectric value if a reaction field cutoff is used,
e.g. |
integrator |
|
friction |
Friction value for the integrator, in inverse time, e.g.
|
platform |
Any valid OpenMM platform string, e.g. |
precision |
Any valid OpenMM platform precision value, e.g.
|
pressure |
Any pressure value, e.g. |
space |
Space in which the simulation should be conducted, e.g. sr.vol.Cartesian |
temperature |
Any temperature value, e.g. |
timestep |
Time between integration steps, e.g.
|
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 }
)