Two interacting point charges.

pme-electrostatics.py

from simtk.openmm import *
from pylab import *

system = System()
system.addParticle(0)
system.addParticle(0)
system.setDefaultPeriodicBoxVectors([10,0,0], [0, 10, 0], [0, 0, 10])
force = NonbondedForce()
force.setNonbondedMethod(NonbondedForce.PME)
force.addParticle(1.0, 1.0, 0.0)
force.addParticle(-1.0, 1.0, 0.0)
force.setCutoffDistance(5.0)
system.addForce(force)

context = Context(system, VerletIntegrator(0))

displacements = linspace(0.5, 2, 1000)
energies = []

for r in displacements:
    context.setPositions([[0,0,0], [r,0,0]])
    energies.append(context.getState(getEnergy=True).getPotentialEnergy()._value)

plot(displacements, energies, 'bx-')
title('PME Interaction Energy. Platform=%s' % context.getPlatform().getName())
xlabel('Displacement')
ylabel('Energy')
show()

woo!