OPLS-WaterSphere-OpenMM.py

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *

pdb = PDBFile('waterSphere.pdb')
#forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
forcefield = ForceField('tip3p.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=NoCutoff)

forces = { system.getForce(index).__class__.__name__ : system.getForce(index) for index in range(system.getNumForces()) }
nonbonded_force = forces['NonbondedForce']
lorentz=CustomNonbondedForce('4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=sqrt(sigma1*sigma2); epsilon=sqrt(epsilon1*epsilon2)')
lorentz.addPerParticleParameter('sigma')
lorentz.addPerParticleParameter('epsilon')
system.addForce(lorentz)


for index in range(nonbonded_force.getNumParticles()):
        charge,sigma,epsilon=nonbonded_force.getParticleParameters(index)
        lorentz.addParticle([sigma,epsilon])
        nonbonded_force.setParticleParameters(index, charge, sigma, epsilon*0)

force = CustomExternalForce('10*max(0, r-1)^2; r=sqrt(x*x+y*y+z*z)')
for i in range(system.getNumParticles()):
    force.addParticle(i, ())
system.addForce(force)
integrator = LangevinIntegrator(1000*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.reporters.append(PDBReporter('output.pdb', 100))
simulation.step(5000)