Pure Liquid Simulations with NAMD and OpenMM

One can use LigParGen server to get parameter and topology files for OpenMM and NAMD MD simulation software. Once you have .pdb .prm .xml and .rtf you can set up you Pure liquid simulations using BOX_MAKER.py program.

python BOX_MAKER.py -p OCT.pdb -b 40

This command creates a 40x40x40 Angstrom box of octanol. This PDB file will then needs to be opened by VMD to create automatic PSF and PDB files for performing MD simulations with NAMD. For OpenMM it can read this file directly and perform minimization and subsequent NPT equilibration. Script used for equlibration and NPT run is given below.

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

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

pdb = PDBFile('box_50_OCT.pdb')
modeller = Modeller(pdb.topology, pdb.positions)
forcefield = ForceField('OCT.xml')
system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer)
system.addForce(MonteCarloBarostat(1*bar,300*kelvin))
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
print('MINIMIZATION STARTED')
simulation = Minimize(simulation,1000)
simulation = Minimize(simulation,1000)
simulation = Minimize(simulation,1000)
simulation = Minimize(simulation,1000)
print('MINIMIZATION DONE')
simulation.reporters.append(PDBReporter('output.pdb', 1000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True,density=True))
simulation.step(100000)
np_equ_pos = simulation.context.getState(getPositions=True).getPositions() 
PDBFile.writeFile(simulation.topology, np_equ_pos,open('NPT_EQ_FINAL.pdb', 'w'))