OpenMM Tutorial

MD simulations with OpenMM

OpenMM is a hardware independent molecular simulation library developed by Pande group at Stanford. OpenMM core libraries are written in C++ but a python wrapper is provided to make the interaction between user and library smoother.

Please look at the below references for more information

• Instaling OpenMM

• JCTC paper on implementation and capabilities

• Documentation

Please install OpenMM using from source code or using conda build if anaconda python is available in your computer. After successful installation continue to the do the tutorial below.

---

Improtant Notes

1. Combination Rules

Its important to note that most of the current force fields except OPLS-AA use Lorentz- Berthelot combination rule $$ \sigma_{ij} = \frac{\sigma_i + \sigma_j}{2} $$ $$\epsilon_{ij} = \sqrt{\epsilon_i \epsilon_j} $$ and is the only combination rule implemented in OpenMM. If one wishes to use geometric combination rule used in OPLS-AA force field, include the function below in your MD codes and call the function to change the combination rule for LJ interactions.

def OPLS_LJ(system):
    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.setNonbondedMethod(nonbonded_force.getNonbondedMethod())
    lorentz.addPerParticleParameter('sigma')
    lorentz.addPerParticleParameter('epsilon')
    lorentz.setCutoffDistance(nonbonded_force.getCutoffDistance())
    system.addForce(lorentz)
    LJset = {}
    for index in range(nonbonded_force.getNumParticles()):
        charge, sigma, epsilon = nonbonded_force.getParticleParameters(index)
        LJset[index] = (sigma, epsilon)
        lorentz.addParticle([sigma, epsilon])
        nonbonded_force.setParticleParameters(
            index, charge, sigma, epsilon * 0)
    for i in range(nonbonded_force.getNumExceptions()):
        (p1, p2, q, sig, eps) = nonbonded_force.getExceptionParameters(i)
        # ALL THE 12,13 and 14 interactions are EXCLUDED FROM CUSTOM NONBONDED
        # FORCE
        lorentz.addExclusion(p1, p2)
        if eps._value != 0.0:
            #print p1,p2,sig,eps
            sig14 = sqrt(LJset[p1][0] * LJset[p2][0])
            eps14 = sqrt(LJset[p1][1] * LJset[p2][1])
            nonbonded_force.setExceptionParameters(i, p1, p2, q, sig14, eps)
    return system

2. Scaling factors for 1-4 interactions

OPLS-AA uses scaling factor of 0.5 for both Lennard Jones and electrostatics, for 1-4 interactions. One should make sure that both solvent and solute have the same 1-4 scale factors. If you find them different, edit the line shown below in forcefield.xml file.

<NonbondedForce coulomb14scale="0.5" lj14scale="0.5">

MD Simulations of Aniline

• Upload the MDL .mol or .pdb file of Aniline or paste SMILES code ('C1=CC=C(C=C1)N') from ChemDraw
• Download the UNK.pdb and UNK.xml files.

Gas-phase minimization

For Gas phase MD simulations, a small script has been written and is shown below. Save it as minimize_gas.py.

For minimization NoCutoffs option is used and no constraints are set on H-bonds.

## Created by Leela S. Dodda for LigParGen Tutorials
## Date Aug 5, 2016
from simtk.openmm import app, KcalPerKJ
import simtk.openmm as mm
from simtk import unit as u
from sys import stdout, exit

def Minimize(simulation,iters=0):
    simulation.minimizeEnergy(maxIterations=iters)
    position = simulation.context.getState(getPositions=True).getPositions()
    energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
    app.PDBFile.writeFile(simulation.topology, position,
                          open('gasmin.pdb', 'w'))
    print 'Energy at Minima is %3.3f kcal/mol' % (energy._value * KcalPerKJ)
    return simulation

temperature = 298.15 * u.kelvin
pdb = app.PDBFile('Aniline.pdb')

modeller = app.Modeller(pdb.topology, pdb.positions)
forcefield = app.ForceField('Aniline.xml')

system = forcefield.createSystem(
    modeller.topology, nonbondedMethod=app.NoCutoff,  constraints=None)
system = OPLS_LJ(system)
integrator = mm.LangevinIntegrator(
    temperature, 1 / u.picosecond,  0.0005 * u.picoseconds)
simulation = app.Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)
simulation = Minimize(simulation,100)

save the above code to minimize_gas.py and run it using below command

bash$ python minimize_gas.py

Gas-phase MD simulation

Once system is minimised, you can add the code below to minimize_gas.py and save it as MD_gas.py to do MD simulations.

simulation.reporters.append(app.PDBReporter('gas_output.pdb', 1000))
simulation.reporters.append(app.StateDataReporter('data.txt', 1000, progress=True, temperature=True, potentialEnergy=True, density=True,totalSteps=10000,speed=True))
simulation.step(100000)