In [1]:
import numpy as np
from simtk import openmm as mm
from simtk.openmm import app
from simtk.unit import *
Set some constants. For this example, we want to set up a system with only coulomb interactions. The OpenMM NonbondedForce computes both the columb and LJ interaction, so to have a system containing only a coulomb force, we can simply set the LJ well depth parameter to zero.
Mathematically, $ V_{LJ} = 4\varepsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] $. Thus if $\varepsilon=0$, then $V_{LJ} = 0$.
In [2]:
N_PARTICLES = 3
MASS = 1.0*dalton
CHARGE = 1.0*elementary_charge
LJ_SIGMA = 1.0*nanometer
LJ_EPSILON = 0.0*kilojoule_per_mole
To get started, create a new System from scratch. The System object is basically a full description of the Hamiltonian of the physical system you'll be simulating. We have to tell the system first how many particles we'll be simulating, and what their masses are.
In [3]:
system = mm.System()
for i in range(N_PARTICLES):
system.addParticle(MASS)
# In case you don't want to look up the documentation, you can get the help
# text associated with each method like this.
print(help(system.addParticle))
Each System can contain zero or more forces, which specify individual terms in the Hamiltonian. There are Forces for the bonded interactions, like HarmonicBondForce, the nonbonded interactions, etc.
We'll add a NonbondedForce, which computes the coulomb and leonard jones interaction.
In [4]:
force = mm.NonbondedForce()
force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
# The force needs to be told about the charge and leonard jones parameters for each of the particles
for i in range(N_PARTICLES):
force.addParticle(CHARGE, LJ_SIGMA, LJ_EPSILON)
# Now that we've created the force, we can add it to our system.
system.addForce(force)
Out[4]:
OpenMM's XMLSerializer can take a System object and save a full description of it to XML. This is pretty useful just so we can check out what's in our system, to make sure our force looks right
In [5]:
print mm.XmlSerializer.serialize(system)
To actually run our simulation, we two more things. The first is an integrator, which specifies some way to evolve the system over time. The second is a Context.
In OpenMM, the Context is a class that handles the actual device that we're running on. It manages the memory, and basically controls the actual physical instantiation of our system and integrator on specific hardware.
In [6]:
integrator = mm.VerletIntegrator(0.002 * picosecond)
context = mm.Context(system, integrator)
print('Running on', context.getPlatform().getName())
The System and Integrator don't actually know about any specific set of positions or velocity. But that is information that the context needs to know. Let's just give it some random positions
In [7]:
context.setPositions(np.random.randn(N_PARTICLES, 3))
context.setVelocities(np.zeros((N_PARTICLES, 3)))
for i in range(10):
integrator.step(10)
state = context.getState(getEnergy=True)
print state.getPotentialEnergy()