OpenMM is a fantastic high performance toolkit for molecular simulation. It is perhaps best known for providing a flexible Python scripting environment for defining and running custom GPU-accelerated molecular dynamics simulations.
In the 2023.2 release of sire we added new sire.convert functions. In the last post we showed how these functions could be used to interconvert with RDKit, thereby adding both smiles and 2D visualisation support.
In this post, we will show how sire.convert can be used to interconvert with OpenMM. This lets sire be used as an interface for setting up OpenMM simulations. We have also added a higher-level set of functions that lets you easily run minimisation and molecular dynamics simulations directly from a set of sire molecules.
For example, here, we will load a molecular system, perform minimisation, and then run molecular dynamics.
import sire as sr
# Load the molecules from our "tutorial" website
mols = sr.load(
sr.expand(sr.tutorial_url, "ala.top", "ala.crd"), silent=True)
# Perform minimisation
mols = mols.minimisation().run().commit()
# Now perform 10 ps of dynamics, saving a snapshot every 0.5 ps
mols = mols.dynamics(
timestep=4*sr.units.femtosecond
).run(10*sr.units.picosecond,
save_frequency=0.5*sr.units.picosecond).commit()
We can view a movie of the resulting trajectory using our integration with NGLView.
mols.view()
We can plot energies of interaction, e.g. between the solute and the water for each of the saved frame of this trajectory using
mols[0].trajectory().energy(mols["water"]).pretty_plot()
This trajectory is, admittedly, very short, on a small molecular system. But the same code would work for larger systems over longer timescales. We just needed something here that you could run easily on our free notebook server ;-)
Under the hood, this works because sire.convert has been called to create an OpenMM system from the loaded sire molecules. This was then manipulated using the standard OpenMM functions to run minimisation and dynamics.
For example, lets create the OpenMM system from our molecules.
omm = sr.convert.to(mols, "openmm")
omm
<openmm.openmm.Context; proxy of <Swig Object of type 'OpenMM::Context *' at 0x7fe17e2f5f90> >
The result is an openmm.Context. This contains a combination of an openmm.System, openmm.Integrator and openmm.Platform.
You can extract these three object using the standard OpenMM functions, e.g.
omm_system = omm.getSystem()
omm_system
<openmm.openmm.System; proxy of <Swig Object of type 'OpenMM::System *' at 0x7fe1b84a6060> >
omm_integrator = omm.getIntegrator()
omm_integrator
<openmm.openmm.VerletIntegrator; proxy of <Swig Object of type 'OpenMM::VerletIntegrator *' at 0x7fe1b84a54b0> >
omm_platform = omm.getPlatform()
omm_platform
<openmm.openmm.Platform; proxy of <Swig Object of type 'OpenMM::Platform *' at 0x7fe17c3b5b40> >
You can combine the omm_system with your own choices of integrator or platform using the standard OpenMM API, if you don’t like the choices made by sire.
To keep the script simple, sire chooses appropriate default settings depending on the molecular system being modelled and the requested simulation timestep. You have complete control of these defaults, and can override any of them using a simple set of options.
The way you do this is by manipulating the Dynamics object which is returned by the .dynamics() function.
d = mols.dynamics(timestep=4*sr.units.femtosecond)
d
Dynamics(completed=0)
You can query the current options by calling the .info(), .constraint() and .ensemble() functions
d.info(), d.constraint(), d.ensemble()
(ForceFieldInfo(
space=PeriodicBox( ( 26.6843, 28.9808, 24.8784 ) ),
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 }
),
'bonds-h-angles',
microcanonical (NVE) ensemble)
In this case, you can see that sire has detected that the molecules are in a periodic box, and so has switched on PME boundary conditions and set a cutoff of 7.5 Å.
It has seen that you haven’t specified a temperature or pressure, so it chooses an integrator appropriate for the microcanonical (NVE) ensemble.
And it uses your choice of a larger timestep (4 fs) to automatically add constraints to all bonds, and to all angles involving hydrogen (the bonds-h-angles constraints).
You can easily change these options, by passing in arguments as described in our detailed guide.
d = mols.dynamics(timestep=1*sr.units.femtosecond,
map={"cutoff": 10*sr.units.angstrom,
"temperature": 25*sr.units.celsius,
"cutoff_type": "REACTION_FIELD"})
d.info(), d.constraint(), d.ensemble()
(ForceFieldInfo(
space=PeriodicBox( ( 26.6843, 28.9808, 24.8784 ) ),
cutoff_type=REACTION_FIELD,
cutoff=10 Å,
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 }
),
'none',
canonical (NVT) ensemble { temperature = 298.15 C })
In this case, we have now switched to a 10 Å reaction field cutoff, and have an integrator that is suited for the canonical (NVT) ensemble at 25°C. Bond and angle constraints have been automatically disabled as a small timestep is being used (you can force a choice by setting the constraint option).
We can then run the simulation by calling d.run() passing in the amount of time we want to simulate, and the amount of time between saving snapshots. Sire then automatically works out the required number of simulation steps based on this, and the desired timestep.
d.run(2*sr.units.picosecond,
save_frequency=0.1*sr.units.picosecond)
Dynamics(completed=6012 ps, energy=-5203.3 kcal mol-1, speed=18.9 ns day-1)
We can call this .run() function as often as we want, changing the amount of time to simulate and the save_frequency as we desire. The coordinates will be appended to our now ever-growing trajectory.
d.run(5*sr.units.picosecond,
save_frequency=1*sr.units.picosecond)
Dynamics(completed=6017 ps, energy=-4410.73 kcal mol-1, speed=20.4 ns day-1)
To finish, we commit this simulation back, converting it back into a sire molecular system. This can be viewed or processed, just as if it had been loaded from a file.
mols = d.commit()
mols.view()
mols[0].trajectory().energy(mols["water"]).pretty_plot()
Notice how the different values of save_frequency for our there trajectories appear differently in the movie and in the above energy graph.
You are not limited to simulating all the molecules. You can simulate any molecule. For example, here we will run MD on the solute (first molecule) loaded from another file.
mols = sr.load(sr.expand(sr.tutorial_url,
"ala.top", "ala.crd"), silent=True)
mol = mols[0].dynamics(
timestep=4*sr.units.femtosecond
).run(
10*sr.units.picosecond,
save_frequency=0.1*sr.units.picosecond
).commit()
mol.view()
mol.trajectory().energy().pretty_plot()
This is just the beginning of what we have planned using our integration with OpenMM. We have better trajectory support which will be part of our 2023.3.0 release later this month, and plans for even more interesting functionality that should be ready for 2023.4.0 later this year. Keep an eye on our GitHub repo if you want to know more ;-)
As with all of our diaries, if you want to try this yourself, please feel free connect to try.openbiosim.org and starting a notebook. You can download the notebook used to generate this post onto the server by running this command in one of the notebook code cells.
! wget https://github.com/OpenBioSim/posts/raw/main/sire/003_openmm/openmm.ipynb
Have a play and let us know what you think. Look forward to our next blog post, where we will show you how to use the new code that will be in 2023.3.0 to load, save, edit and convert molecular trajectories.
