Problem with molecule going out of water box

[Jonathan Saboury]

After about a 720 picosecond simulation of a prmtop and inpcrd file created by AMBER from a pdb created in spartan '08 in OpenMM I get a pdb file showing the simulation that is visualized with VMD.

At 360 picoseconds the molecule seems to go out of the box of water. I was told before that this was an imaging issue and everything is fine, but I have some doubts of this. When the molecule is inside water there are strong pi-pi interactions between benzenes, but as soon as it goes out of the water box those interactions are destroyed.

Here is the commands to get the prmtop and inpcrd files: http://pastebin.com/raw.php?i=63sFFQ6q

Thank you for your time, it is very much appreciated!

[Peter Eastman]

Hi Jonathan,

What nonbonded method are you specifying when you call createSystem()? (Can you post your script?)

Does the PDB file generated by OpenMM contain a CRYST1 record?

Peter

[Jonathan Saboury]

The .py I used is:

Code: Select all

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


prmtop = AmberPrmtopFile('2cNitro.prmtop')
inpcrd = AmberInpcrdFile('2cNitro.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('outputPyMol.pdb', 1000))

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


simulation.reporters.append( StateDataReporter('outputPyMolEnergy.txt', 1000, step=True, potentialEnergy=True, temperature=True) )

simulation.step(750000)

And yes there is, it is at the start of the text file. For this case I used a waterbox that is 30 Angstroms, which let me simulate the molecule longer getting more interesting results at the cost of much more computations

I don't know if there are more after it because it is taking a long time to search (3 million lines):

Code: Select all

CRYST1   85.960   76.851   98.363  90.00  90.00  90.00 P 1           1 
MODEL        1
ATOM      1  C   UNK A   1      47.802  40.151  45.490  1.00  0.00
ATOM      2  C1  UNK A   1      47.562  40.203  48.250  1.00  0.00

The AMBER people are telling me that the box itself should change coordinates/move around as well, is this how the box works in OpenMM as well?

Thanks!

[Peter Eastman]

Ok, periodic boundary conditions are definitely being applied. Periodic boundary conditions don't affect the position of any atom, so there's nothing to prevent atoms from drifting beyond what you think of as the "box". They affect how forces are calculated, so when it computes the vector between two atoms, it's computed modulo the box dimensions.

I gather you're using OpenMM 4.1? Try upgrading to 5.0. In 5.0, the PDB writer also applies periodic boundary conditions, so all your molecules get translated into a single periodic box. This has nothing to do with how the simulation is run, but it makes your output look more like you expect.

Peter

[Jonathan Saboury]

Yes was running 4.1, will be updating to 5.0 to see how it changes the output, thank you!

[Jonathan Saboury]

Sorry for the late response. 5.0 did not change the molecule from going out of periodic box.

Perhaps I did not install 5.0 correctly? But doubt it, I deleted the old 4.1 installation folder and put in 5.0.