Minimizing and running a simulation of an organic molecule

[Peter Eastman]

Exactly. PDBFile will load your file whatever is in it, but ForceField won't know how to create a System from it if there's a molecule that it doesn't know about.

This assumes, of course, that you know what force field parameters you want to use for your molecules. If not, an alternative is to use AmberTools to generate a prmtop file, since it can generate force field parameters for arbitrary small molecules.

Peter

[Jonathan Saboury]

Installed Ubuntu and AmberTools12. I have started XLEaP, loaded the .pdb I had, but when I try to generate the .prmtop and .inpcrd it gives an error.

I know that this question belongs to the AMBER mailing list. I have posted a couple times there, but I only see my first question which does not apply to me anymore. I have a general idea of how mailing lists work but not 100% sure my question is being sent through (as it isn't being posted in the archives like the first one).

Hopefully you will know how to deal with this problem, but if not, it is fine. You went beyond the call of duty helping me and it is very much appreciated.

As a side note, once I have figured out how to do this i'll post my solutions because undoubtedly someone else may have the same problem and have been dissuaded due to difficulty.

Here is my problem, I am so close I can taste it!

Thanks again!

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

I loaded this .pdb (http://pastebin.com/9ZMS3E6E) with "hexane1 = loadpdb "/home/jonathan/amber12/PDBs/hexane.pdb"

and it outputs:

Code: Select all

Loading PDB file: /home/jonathan/amber12/PDBs/hexane.pdb
Unknown residue: PSD    number: 0    type:   Terminal/last
..relaxing end constraints to try for a dbase match
  -no luck
Creating a new UNIT for residue: PSD sequence:   1
Created a new atom named: C01 within residue:   .R<PSD 1>
...
Created a new atom named: H14 within residue:   .R<PSD 1>
    total atoms in file: 20
    The file contained 20 atoms not in residue templates

Then I type the command "saveamberparm hexane1 hexane.prmtop hexane.inpcrd"
and it outputs:

Code: Select all

Checking Unit.
FATAL"    Atom .R<PSD 1>.A<C01 1> does not have a type.
....
FATAL"    Atom .R<PSD 1>.A<H14 20> does not have a type.
Failed to generate parameters
Parameter file was not saved.

How do I make it save the prmtop and inpcrd for this file? Thanks!

[Jonathan Saboury]

Also if I entered "edit hexane1" it would visualize the .pdb correctly. So I know it read it correctly.

[Morgan Lawrenz]

Hello,
Because you are working with a non-protein residue, here you need to parametrize your molecule using antechamber and GAFF, which should be pretty simple for a hexane molecule.
The tutorial here gives a step-by-step of the files you need to generate and load into Leap:
http://ambermd.org/tutorials/basic/tutorial4b/
Thanks

[Jonathan Saboury]

Thank you Morgan,

I have done as the tutorial instructed, created the inpcrd and prmtop files (which is in the zip attached) and executed the below code:

Code: Select all

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


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


system = prmtop.createSystem(nonbondedMethod=PME,nonbondedCutoff=1*nanometer, constraints=HBonds)

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', 1000))

simulation.reporters.append( StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True) )

simulation.step(10000)

But I get this error:

Code: Select all

C:\Program Files (x86)\OpenMM\hexane>python hexane.py
Traceback (most recent call last):
  File "hexane.py", line 11, in <module>
    system = prmtop.createSystem(nonbondedMethod=PME,nonbondedCutoff=1*nanometer
, constraints=HBonds)
  File "C:\Python27-32Bit\lib\site-packages\simtk\openmm\app\amberprmtopfile.py"
, line 137, in createSystem
    raise ValueError('Illegal nonbonded method for a non-periodic system')
ValueError: Illegal nonbonded method for a non-periodic system

Thanks again!

[Peter Eastman]

This means that when you set up your system in AmberTools, you didn't tell it to use periodic boundary conditions.

Peter

[Jonathan Saboury]

Thank you Peter,

I set periodic boundary conditions with Amber12 XLEaP.

Was having problems with setting the right "nonbondedCutoff". If i set it to "0.5xnanometer" it was larger than half the system. Set it too small (0.01 x nanometer) got a gpu related error and I suspect it was just calling for too much memory. Set it to "0.1*nanometer" and presto, fixed that error but got another one.

I have said all of this because I'm curious what it means. Does it mean that that atom is only affected by other atoms that are 0.1*nanometer's away?

Anyways, OpenMM is having trouble outputting to the .pdb, here is the error I am getting currently, and attaching the input files.

Code:

Code: Select all

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


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


system = prmtop.createSystem(nonbondedMethod=PME,nonbondedCutoff=0.1*nanometer, constraints=HBonds)

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', 1000))

simulation.reporters.append( StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True) )

simulation.step(10000)

Error:

Code: Select all

c:\hexane>python hexane.py
Traceback (most recent call last):
  File "hexane.py", line 25, in <module>
    simulation.step(10000)
  File "C:\Python27-32Bit\lib\site-packages\simtk\openmm\app\simulation.py", lin
e 127, in step
    reporter.report(self, state)
  File "C:\Python27-32Bit\lib\site-packages\simtk\openmm\app\pdbreporter.py", li
ne 78, in report
    PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, sel
f._nextModel)
  File "C:\Python27-32Bit\lib\site-packages\simtk\openmm\app\pdbfile.py", line 2
45, in writeModel
    raise ValueError('Particle position is NaN')
ValueError: Particle position is NaN

As always, thank you for your time.

[Peter Eastman]

Hi Jonathan,

Are you sure you set your periodic box size correctly? 0.1 nanometer is roughly the width of a single hydrogen atom. A typical nonbonded cutoff is around 1 nanometer. What is your box size supposed to be?

Peter

[Jonathan Saboury]

Well, I have been using the LEaP command "setBox HEX centers" (thinking it will do what you ask automatically) where HEX is from the mol2 file created form antechamber.

I honestly have no idea what most of the commands do so I think the best thing I can do for you to grasp what I am doing wrong is post the command prompt i/o.

I put "hexane.pdb" created from PyMol into the directory "/home/jonathan/".

I first put "export AMBERHOME=/home/jonathan/amber12" which outputs nothing.

then "export AMBERHOME=/home/jonathan/amber12" which outputs:

Code: Select all

Running: /home/jonathan/amber12/bin/bondtype -j full -i ANTECHAMBER_BOND_TYPE.AC0 -o ANTECHAMBER_BOND_TYPE.AC -f ac

Running: /home/jonathan/amber12/bin/atomtype -i ANTECHAMBER_AC.AC0 -o ANTECHAMBER_AC.AC -p gaff
Total number of electrons: 50; net charge: 0

Running: /home/jonathan/amber12/bin/sqm -O -i sqm.in -o sqm.out

Running: /home/jonathan/amber12/bin/am1bcc -i ANTECHAMBER_AM1BCC_PRE.AC -o ANTECHAMBER_AM1BCC.AC -f ac -p /home/jonathan/amber12/dat/antechamber/BCCPARM.DAT -s 2 -j 1

then I run "$AMBERHOME/bin/parmchk -i hexane.mol2 -f mol2 -o hexane.frcmod" which outputs nothing.

then start tleap with "$AMBERHOME/bin/tleap -f leaprc.ff99SB"
outputs this:

Code: Select all

-I: Adding /home/jonathan/amber12/dat/leap/prep to search path.
-I: Adding /home/jonathan/amber12/dat/leap/lib to search path.
-I: Adding /home/jonathan/amber12/dat/leap/parm to search path.
-I: Adding /home/jonathan/amber12/dat/leap/cmd to search path.
-f: Source leaprc.ff99SB.

Welcome to LEaP!
(no leaprc in search path)
Sourcing: /home/jonathan/amber12/dat/leap/cmd/leaprc.ff99SB
Log file: ./leap.log
Loading parameters: /home/jonathan/amber12/dat/leap/parm/parm99.dat
Reading title:
PARM99 for DNA,RNA,AA, organic molecules, TIP3P wat. Polariz.& LP incl.02/04/99
Loading parameters: /home/jonathan/amber12/dat/leap/parm/frcmod.ff99SB
Reading force field modification type file (frcmod)
Reading title:
Modification/update of parm99.dat (Hornak & Simmerling)
Loading library: /home/jonathan/amber12/dat/leap/lib/all_nucleic94.lib
Loading library: /home/jonathan/amber12/dat/leap/lib/all_amino94.lib
Loading library: /home/jonathan/amber12/dat/leap/lib/all_aminoct94.lib
Loading library: /home/jonathan/amber12/dat/leap/lib/all_aminont94.lib
Loading library: /home/jonathan/amber12/dat/leap/lib/ions94.lib
Loading library: /home/jonathan/amber12/dat/leap/lib/solvents.lib

then run "source leaprc.gaff" which outputs:

Code: Select all

----- Source: /home/jonathan/amber12/dat/leap/cmd/leaprc.gaff
----- Source of /home/jonathan/amber12/dat/leap/cmd/leaprc.gaff done
Log file: ./leap.log
Loading parameters: /home/jonathan/amber12/dat/leap/parm/gaff.dat
Reading title:
AMBER General Force Field for organic molecules (Version 1.4, March 2010) add. info. at the end

then I load the mol2 created with "HEX = loadmol2 hexane.mol2" which outputs:

Code: Select all

Loading Mol2 file: ./hexane.mol2
Reading MOLECULE named PSD

then I run "check HEX" which outputs:

Code: Select all

Checking 'HEX'....
Checking parameters for unit 'HEX'.
Checking for bond parameters.
Checking for angle parameters.
Unit is OK.

then i run "loadamberparams hexane.frcmod" which outputs:

Code: Select all

Loading parameters: ./hexane.frcmod
Reading force field modification type file (frcmod)
Reading title:
remark goes here

then i run "saveoff HEX hexane.lib" which outputs:

Code: Select all

 Creating hexane.lib
Building topology.
Building atom parameters.

then i run "setBox HEX centers" which outputs:

Code: Select all

Box dimensions:  6.542000 4.919000 3.802000

then I run "saveamberparm HEX hexane.prmtop hexane.inpcrd" which outputs:

Code: Select all

Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
 total 0 improper torsions applied
Building H-Bond parameters.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
  (Residues lacking connect0/connect1 - 
   these don't have chain types marked:

	res	total affected

	PSD	1
  )
 (no restraints)

Then I load those two files in the .py in the previous post.

I have posted all the files created during this process in case your need them.

As always, thank you so much for your time.

[Jonathan Saboury]

Figiured out I needed to use "solvateBox" instead.

Was able to run a simulation of Hexane

Thank you! (not completely done yet, having a problem with compounds with carbonyls, but this is a AMBER issue)