OpenMM5.1 'first periodic box must be parallel to x' problem

[Krishna dev Oruganty]

Hello,

I am trying to simulate a protein (based on 1mq4, 272 residues) with OpenMM 5.1 compiled from source (compilation gave no errors, and all CUDA tests passed). Hardware specs are dual X5560 CPU and a GTX670 with system memory of 24GB. My system is running Ubuntu 12.04 and I have python 2.7.3, CUDA 5.0 with NVIDIA driver 310.19, and GCC v4.6.3. I used the builder script from your site (https://openmm.herokuapp.com/) to start the simulation (implicit solvent, AMBER99sb-ildn forcefield). I had modeled the missing residues in the PDB using Modeller package (salilab) and added hydrogens using OpenMM modeller class.

At random times, the simulation crashes. The first time, it happened after 10,000 steps and later after 36,000 steps and the final one at 6000 steps. I have simulated the same system with explicit water too, but the results are the same (with random crashes). I tried openCL also, but it keeps crashing. I also tried CUDA single precision and double precision, but the problem persists. There are two things I want to ask:

1) If the simulation is setup for implicit solvent, why do we need a periodic box (see error message below) ? What am I doing wrong here ?
2) Is this a problem with the GPU returning garbage values and that crashes the simulation ? Apart from dumping every frame into a PDB file, is there any other way to check if this is happening ? And how can I improve the stability of the simulations ?

The output when running the simulation is as follows :

Minimizing...
Equilibrating...
Running Production...
#"Step","Potential Energy (kJ/mole)","Temperature (K)"
1000,-30211.8785067,283.386134539
2000,-29293.3160687,295.469092923
3000,-28943.8682014,305.244706957
4000,-28931.4746189,306.729539459
5000,-29409.7607804,300.015820269
6000,-29220.6446371,296.14181751
Traceback (most recent call last):
File "openmm.py", line 42, in <module>
simulation.step(100000)
File "/usr/local/lib/python2.7/dist-packages/simtk/openmm/app/simulation.py", line 105, in step
self.integrator.step(10) # Only take 10 steps at a time, to give Python more chances to respond to a control-c.
File "/usr/local/lib/python2.7/dist-packages/simtk/openmm/openmm.py", line 9433, in step
return _openmm.LangevinIntegrator_step(self, *args)
Exception: First periodic box vector must be parallel to x.

I have used the same starting structure with Gromacs 4.6.1 and AMBER99sb-ILDN+TIP3P water and I am able to get stable simulation till 200ns. So, I dont think the structure is the problem.

The script I am using is as follows:

Code: Select all

##########################################################################
# this script was generated by openmm-builder. to customize it further,
# you can save the file to disk and edit it with your favorite editor.
##########################################################################

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

pdb = PDBFile('Aur_1mq4.pdb')
forcefield = ForceField('amber99sbildn.xml', 'amber99_obc.xml')

modeller = Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield)


system = forcefield.createSystem(modeller.topology, nonbondedMethod=NoCutoff,
     constraints=HBonds, rigidWater=False, removeCMMotion=True)
integrator = LangevinIntegrator(300*kelvin, 1.0/picoseconds, 2.0*femtoseconds)
integrator.setConstraintTolerance(0.00001)
system.addForce(MonteCarloBarostat(1*atmospheres, 300*kelvin, 25))

platform = Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed'}
simulation = Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)

print('Minimizing...')
simulation.minimizeEnergy()

simulation.context.setVelocitiesToTemperature(300*kelvin)
print('Equilibrating...')
simulation.step(100)

simulation.reporters.append(PDBReporter('output.pdb', 1000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True,
    potentialEnergy=True, temperature=True))

print('Running Production...')
simulation.step(100000)
print('Done!')

[Peter Eastman]

Hi Krishna,

It sounds to me like your simulation is blowing up, probably from using settings that are just outside of (or right on the edge of) the stable region. Looking at your script I notice this:

system = forcefield.createSystem(modeller.topology, nonbondedMethod=NoCutoff,
constraints=HBonds, rigidWater=False, removeCMMotion=True)
integrator = LangevinIntegrator(300*kelvin, 1.0/picoseconds, 2.0*femtoseconds)

For explicit solvent, that definitely won't be stable. With flexible water (you specified rigidWater=False), you need a smaller time step than 2 fs. 0.5 fs is a typical choice in that case. But do you really want flexible water? Most water models are parametrized for rigid water.

If you don't have any explicit water then the settings are more reasonable, but they may still be going a little too close to the edge. Try monitoring the temperature during your simulation. It will fluctuate around some average, but that average should be very close to 300K. If it's more than two or three degrees above that, you really should either decrease your time step slightly or increase your friction.

If you do that, does the error go away?

Peter

[Krishna dev Oruganty]

Hello Peter,

I tried with time steps of 0.5fs, 1.0fs and 2.0fs with friction terms of 1.0/ps and 2.0/ps. But it kept crashing randomly. But, if I increase the equilibration from 100 steps to 1000 steps, then the simulation completes 100000 steps of production run without problem even with the original parameters (2.0 fs time step and 1.0/ps friction term). The commands for equilibration and production runs are the same, so I had not given that much importance earlier. But, I guess adding the PDBReporter statement may have changed something.

If we have implicit solvent, then the program should ignore the rigid water directive, isn't it ? Also, if it is an implicit solvent simulation, do we need a periodic box ?

-Krishnadev

[Peter Eastman]

Hi Krishnadev,

I realized what the problem is. You're adding a MonteCarloBarostat to the system, which doesn't make sense for a nonperiodic system. Let me explain what's going on.

The Context always defines some set of periodic box vectors. It's up to the individual forces to decide whether to use them or not. Since you specified nonbondedMethod=NoCutoff, they just ignore them, but they're still there.

The barostat keeps adjusting the size of the periodic box, and since you have finite pressure, it mostly tends to make the box smaller. It keeps shrinking the box until the resistance from compressing the system balances out the pressure. But since all the forces just ignore the box size, it never feels resistance. It just keeps shrinking the box until the size reaches zero, at which point you get an error message. (What it actually means by, "First periodic box vector must be parallel to x," is that (0, 0, 0) isn't parallel to (1, 0, 0). Admittedly, in this case that doesn't make it very clear what the root problem is.)

So the short answer: remove the barostat!

Peter

[Krishna dev Oruganty]

Hi Peter,

Thanks a lot for the help. I had not realized that Barostat was still on even in the implicit solvent simulations. Removing it resolved the problem. Perhaps, the script builder can be modified to warn about this scenario.

Thanks,
Krishnadev

[Robert McGibbon]

I've updated the script builder (see https://github.com/rmcgibbo/openmm-webbuilder/issues/8) to issue a warning in this scenario.

[Krishna dev Oruganty]

Hello Robert and Peter,

The script builder works fine now and warns about the problem.

Thanks for all your help.

-Krishnadev