simtk.org

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)

Hi Brian,

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

Thanks,

- Lee-Ping

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

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

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