simtk.org

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

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

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)

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

Hi Ian,

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

Peter

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

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

Peter

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.