LigParGen

OPLS/CM1A Parameter Generator for Organic Ligands

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

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.

Important Notes

1. Combination Rules

Most of the current standard force fields except OPLS-AA use the Lorentz- Berthelot combination rule to estimate epsilon (ε_ij) and sigma (σ_ij) Van Der Waals parameters σ_ij = (σ_i + σ_j)/2 ; ε_ij = √ ε_i ε_j and is the only combination rule implemented in OpenMM. Users must include and call the function below to use the OPLS-AA geometric combination rule in the MD simulation to change the 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 the 0.5 scaling factor for 1-4 interaction terms (Lennard Jones and electrostatics). Please, make sure that solvent and solute have the same 1-4 scale factors in the force field parameter files. If you find them different, edit the line shown below in forcefield.xml file.

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

OpenMM Ligand Simulations

1. Preparing Ligand System

LigParGen server provides OPLS-AAM templates with CM1A/CM1A-LBCC charges for small organic molecules. In general, molecular dynamics simulations are focused on protein/NA-ligand interactions rather than just small molecules but OPLS-AA is not yet implemented in OpenMM templates. For this reason, in this tutorial, just a protocol to perform ligand simulations using LigParGen server will be explained in detail.

OpenMM and python must be installed in your computer to perform this tutorial.

Aniline

Ligand template generation will be performed in the LigParGen server following the next two steps:

Upload a .mol or .pdb file of Aniline or paste SMILES code C1=CC=C(C=C1)N
OPTIONAL: Aniline molecule can be drawn directly in the Draw Molecule section in the LigParGen server
Download the .pdb and .xml files.

2. Gas-phase minimization

First, the system must be minimized to perform a molecular dynamics simulation. It is required even if the structure optimization has been performed in the LigParGen server in order to acomodate the water molecules. In this case, the optimization is optional (if it is done in the server) because a simulation in gas phase will be done.

A small script to perform the minimization is shown below. Copy and paste the code lines in a file called 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)

Run the python code using below command

bash$ python minimize_gas.py

3. Gas-phase MD simulation

Once system is minimized, add the code below to minimize_gas.py and save it as MD_gas.py. Then, to run an MD simulations use the same previous minimization command.

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)