Water collapse in simulation

The functionality of OpenMM will (eventually) include everything that one would need to run modern molecular simulation.
POST REPLY
User avatar
Brian Geiss
Posts: 12
Joined: Thu Dec 30, 2010 8:19 am

Water collapse in simulation

Post by Brian Geiss » Fri Nov 01, 2013 2:43 pm

I recently ran a 50 ns explicit water simulation and had the waters and ions collapse into a tight ball in the center of the box for some unknown reason. Some waters remained associated with my protein, but <99% dropped into this cluster. At the end of the simulation the water ball started expanding, again for unknown reasons. While interesting to watch, it's not that useful to me. I haven't been able to find anything about this kind of behavior and wanted to see if anyone else has had similar issues. I've copied my prep and run scripts below in case I've made a simple mistake in the set up.

Thanks for help in advance.

-Brian

Prep script

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
print('Loading...')
pdb = PDBFile('3EVG_M219Q.pdb')
forcefield = ForceField('amber99sbildn.xml', 'tip3p.xml')
modeller = Modeller(pdb.topology, pdb.positions)
print('Adding hydrogens...')
modeller.addHydrogens(forcefield)
print('Adding solvent...')
modeller.addSolvent(forcefield, model='tip3p', padding=1*nanometer, ionicStrength=0.1*molar)
print('Minimizing...')
system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME, nonbondedCutoff=1.2*nanometer)
integrator = VerletIntegrator(0.002*picoseconds)
simulation = Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)
simulation.minimizeEnergy()
print('Saving...')
positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(simulation.topology, positions, open('3EVG_M219Q_fixed_ildn_TIP3p.pdb', 'w'))
print('Done')


Run script
##########################################################################
# 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
import time;

startticks = time.time()
pdb = PDBFile('3EVG_M219Q_fixed_ildn_TIP3p.pdb')
forcefield = ForceField('amber99sbildn.xml', 'tip3p.xml')

system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1.2*nanometer,
constraints=HBonds, rigidWater=True)

system.addForce(MonteCarloBarostat(1*bar, 310*kelvin))

integrator = LangevinIntegrator(310*kelvin, 1.0/picoseconds, 2.0*femtoseconds)
integrator.setConstraintTolerance(0.00001)

platform = Platform.getPlatformByName('OpenCL')

properties = {'OpenCLPrecision': 'mixed'}
simulation = Simulation(pdb.topology, system, integrator, platform, properties)
simulation.context.setPositions(pdb.positions)

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

simulation.context.setVelocitiesToTemperature(310*kelvin)
print('Equilibrating...')
simulation.step(100000)

simulation.reporters.append(PDBReporter('3EVG_M219Q_50ns_mixed_ildn_tip3p_OpenCL_50psoutput.pdb', 50000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))
simulation.reporters.append(StateDataReporter('3EVG_M219Q_50ns_mixed_ildn_tip3_OpenCL_50psoutput', 50000, step=True,
time=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True,
temperature=True, volume=True, density=True))

print('Running Production...')
simulation.step(50000000)
print('Done!')
endticks = time.time()
final = (endticks-startticks)/60
print ('total run time(minutes)', final)

User avatar
Lee-Ping Wang
Posts: 102
Joined: Sun Jun 19, 2011 5:14 pm

Re: Water collapse in simulation

Post by Lee-Ping Wang » Fri Nov 01, 2013 9:47 pm

Hi Brian,

Your scripts look fine to me. How tight is the ball of solvent? Can you post a screenshot?

Thanks,

- Lee-Ping

User avatar
Brian Geiss
Posts: 12
Joined: Thu Dec 30, 2010 8:19 am

Re: Water collapse in simulation

Post by Brian Geiss » Mon Nov 04, 2013 2:09 pm

Here's a screenshot of Frame 1 and frame 100. The water ball re-expands at the end of the simulation.
Frame 1
Frame1 water.png
Frame 100
Frame 100 water.png

User avatar
Lee-Ping Wang
Posts: 102
Joined: Sun Jun 19, 2011 5:14 pm

Re: Water collapse in simulation

Post by Lee-Ping Wang » Tue Nov 05, 2013 7:26 am

Hi Brian,

If you look near the middle of the ball, are the molecules overlapping (more precisely, do any intermolecular distances drop below 1.4 angstrom or so)? I'm curious about whether certain interaction terms are missing which would cause unphysical geometries to appear.

Also, if you could print out the energy and temperature as this simulation is progressing, that would be ideal.

Thanks,

- Lee-Ping

User avatar
Brian Geiss
Posts: 12
Joined: Thu Dec 30, 2010 8:19 am

Re: Water collapse in simulation

Post by Brian Geiss » Wed Nov 06, 2013 12:07 pm

Hi Lee-Ping,

Thanks for taking the time to look at my data. I've looked attached total energy and temperature graphs below, and there doesn't look like much is changing there. I've also attached a screenshot in the middle of the water ball and have set the H-bond distance to 1.4A. As you can see there are a lot of waters that are within that distance, and if I drop the H-bond distance to 0.5A I can still see a few in the center of the mass. Using the standard 3.0A there are too many H-bonds between waters for VMD to model. Could this be an incompatibility between the TIP3P water model and the amber99sbildn forcefield?

Thanks again,

-Brian
Energy.pdf
(49.43 KiB) Downloaded 76 times
Temperature.pdf
(52.84 KiB) Downloaded 74 times
Water Sphere 1.4A bond.png

POST REPLY