Getting Started with Molecular Dynamics

[Ian Dunn]

I am having some difficulty getting started with OpenMM. My goal is to be able to take any arbitrary organic molecule surrounded by explicit water molecules and to run classical molecular dynamics. I want to generate an animation of the molecular dynamics that can be viewed in jmol. Right now I cannot even properly load my system into OpenMM. Below is the .pdb file of my system that was generated from a .xyz file using Babel and then solvated with explicit waters in Packmol. The error message I get in OpenMM is the following. Does anyone have any thoughts on how I can fix this, or if there is another force field I can use for this problem that will work properly? Thanks!

>>> pdb = PDBFile('solute.pdb')
>>> forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
>>> system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/n/sw/Python-2.7.3/lib/python2.7/site-packages/simtk/openmm/app/forcefield.py", line 314, in createSystem
raise ValueError('No template found for residue %d (%s). This might mean your input topology is missing some atoms or bonds, or possibly that you are using the wrong force field.' % (res.index+1, res.name))
ValueError: No template found for residue 1 (LIG). This might mean your input topology is missing some atoms or bonds, or possibly that you are using the wrong force field.


OpenMM Input:

HEADER
TITLE Built with Packmol
REMARK Packmol generated pdb file
REMARK Home-Page: http://www.ime.unicamp.br/~martinez/packmol
REMARK
HETATM 1 H HOH A 1 -5.987 0.206 1.113
HETATM 2 H HOH A 1 -7.384 1.014 0.866
HETATM 3 O HOH A 1 -6.421 0.952 0.607
HETATM 4 H HOH A 2 -5.631 5.010 1.055
HETATM 5 H HOH A 2 -5.168 6.106 -0.064
HETATM 6 O HOH A 2 -4.848 5.410 0.578
HETATM 7 H HOH A 3 1.091 -3.482 -6.650
HETATM 8 H HOH A 3 0.232 -4.367 -5.580
HETATM 9 O HOH A 3 0.798 -4.405 -6.403
HETATM 10 H HOH A 4 -5.588 -0.405 -5.593
HETATM 11 H HOH A 4 -4.682 0.439 -4.529
HETATM 12 O HOH A 4 -5.485 0.472 -5.123
HETATM 13 H HOH A 5 4.812 2.847 -2.124
HETATM 14 H HOH A 5 5.125 4.263 -2.875
HETATM 15 O HOH A 5 4.722 3.354 -2.981
HETATM 16 H HOH A 6 -4.102 4.802 -2.672
HETATM 17 H HOH A 6 -3.464 6.157 -2.022
HETATM 18 O HOH A 6 -3.666 5.192 -1.861
HETATM 19 H HOH A 7 -6.581 -2.717 -1.716
HETATM 20 H HOH A 7 -6.272 -4.215 -2.289
HETATM 21 O HOH A 7 -6.883 -3.671 -1.714
HETATM 22 H HOH A 8 -3.224 -2.091 -2.829
HETATM 23 H HOH A 8 -3.324 -3.213 -4.012
HETATM 24 O HOH A 8 -3.695 -2.349 -3.673
HETATM 25 H HOH A 9 -6.162 -0.045 -2.079
HETATM 26 H HOH A 9 -7.735 0.361 -1.918
HETATM 27 O HOH A 9 -6.935 -0.007 -1.446
HETATM 28 H HOH A 10 -4.986 3.127 -1.232
HETATM 29 H HOH A 10 -4.721 1.536 -1.488
HETATM 30 O HOH A 10 -5.415 2.224 -1.279
ATOM 31 O LIG 1 2.210 -0.107 -0.400 1.00 0.00 O
ATOM 32 C LIG 1 1.591 0.892 0.459 1.00 0.00 C
ATOM 33 C LIG 1 1.238 0.287 1.826 1.00 0.00 C
ATOM 34 O LIG 1 0.090 -0.550 1.699 1.00 0.00 O
ATOM 35 P LIG 1 0.338 -2.164 1.603 1.00 0.00 P
ATOM 36 O LIG 1 0.964 -2.528 2.940 1.00 0.00 O1-
ATOM 37 O LIG 1 1.277 -2.290 0.410 1.00 0.00 O
ATOM 38 O LIG 1 -1.052 -2.739 1.374 1.00 0.00 O1-
ATOM 39 C LIG 1 2.607 2.066 0.470 1.00 0.00 C
ATOM 40 O LIG 1 3.010 2.253 -0.904 1.00 0.00 O
ATOM 41 C LIG 1 2.093 3.407 1.015 1.00 0.00 C
ATOM 42 O LIG 1 1.710 3.283 2.400 1.00 0.00 O
ATOM 43 C LIG 1 3.168 4.519 0.954 1.00 0.00 C
ATOM 44 O LIG 1 3.669 4.906 2.023 1.00 0.00 O
ATOM 45 C LIG 1 3.529 5.086 -0.425 1.00 0.00 C
ATOM 46 O LIG 1 3.771 6.476 -0.389 1.00 0.00 O
ATOM 47 P LIG 1 4.488 7.090 -1.732 1.00 0.00 P
ATOM 48 O LIG 1 3.700 6.499 -2.886 1.00 0.00 O1-
ATOM 49 O LIG 1 5.930 6.612 -1.639 1.00 0.00 O
ATOM 50 O LIG 1 4.314 8.599 -1.580 1.00 0.00 O1-
ATOM 51 H LIG 1 1.813 -0.994 -0.126 1.00 0.00 H
ATOM 52 H LIG 1 0.679 1.195 -0.067 1.00 0.00 H
ATOM 53 H LIG 1 0.962 1.056 2.547 1.00 0.00 H
ATOM 54 H LIG 1 2.083 -0.276 2.235 1.00 0.00 H
ATOM 55 H LIG 1 3.508 1.757 1.015 1.00 0.00 H
ATOM 56 H LIG 1 2.987 1.339 -1.270 1.00 0.00 H
ATOM 57 H LIG 1 1.213 3.757 0.464 1.00 0.00 H
ATOM 58 H LIG 1 2.474 3.645 2.899 1.00 0.00 H
ATOM 59 H LIG 1 2.683 4.907 -1.096 1.00 0.00 H
ATOM 60 H LIG 1 4.407 4.546 -0.785 1.00 0.00 H
END

[Peter Eastman]

Hi Ian,

Right now, OpenMM doesn't include a tool to generate force field parameters for arbitrary organic molecules. The built in force fields just provide parameters for the standard amino acids, nucleotides, etc. So I suggest using AmberTools to prepare your system. You can use antechamber to generate force field parameters, then save prmtop and inpcrd files that OpenMM can read.

Peter

[Ian Dunn]

Ok, I used antechamber to generate the prmtop and inpcrd files and MD works great with one molecule. How do I then add explicit waters? In particular I want to add exactly 10 explicit waters to the system. My python script is included down below. Since I am not using the Modeller class I do not know how to use the addSolvent method in this situation. Thanks for your help!

Python Script:

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

prmtop = AmberPrmtopFile('solute.prmtop')
inpcrd = AmberInpcrdFile('solute.inpcrd')

system = prmtop.createSystem(nonbondedMethod=NoCutoff, constraints=HBonds)
system.addSolvent(ForceField('amoeba2009.xml'))
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(prmtop.topology, system, integrator)
simulation.context.setPositions(inpcrd.positions)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter('output.pdb', 10))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))
simulation.step(10000)

[Ian Dunn]

Also, please ignore the system.addSolvent line in the script above. When I tried to add that line the script failed.

[Peter Eastman]

Hi Ian,

The simplest thing will be to just build the whole model with AmberTools, water included.

Peter

[Ian Dunn]

Hi Peter,

Thanks for your help. Is there a straightforward way to do this with exactly 10 explicit waters. Also I have never used Amber before so if you get a chance would you mind posting the code I'd need to run to set up this system in AmberTools and run it in OpenMM? I get my initial coordinates in the form of an SDF file and from there I want to be able to solvate the system and run MD on it. Once again I appreciate your help.

Regards,
Ian

[Peter Eastman]

I don't have much experience with AmberTools myself. Anyone with more experience want to comment?

Peter

[Morgan Lawrenz]

Here are some links to a few good Amber Tutorials on setting up simulations with AmberTools.

From Ross Walker (Amber developer):
http://ambermd.org/tutorials/
I recommend looking at both B1 (for simple MD set up) and B4 (for simulating a compound parametrized with GAFF).

From the Rizzo lab, a good summary but uses implicit solvent:
http://ringo.ams.sunysb.edu/index.php/2 ... reptavidin

In all of the above, you generate parameters for the molecule, and can solvate with a pre-equilibrated water box using solvateBox.

Adding just a few waters might be trickier, here's a help forum answer that recommends a specific Amber took ParmEd to strip out extra waters:
http://archive.ambermd.org/201205/0053.html

If you are having Amber-specific problems, I'd recommend posting on their help list, they are also very responsive.
Good Luck.