Water molecules go out of bounding box while running MD

[JSLJ23]

Hi openmm developers,

I was trying to run a simple MD simulation to fold a small peptide of 10 amino acids in length, and loading the PDB file directly with openmm.app.PDBFile but for some reason, after running the simulation for 10,000,000 steps, all the water molecules drift far away from the protein.

Here is what I tried to run:

Imports:

from openmm.app import PDBFile, Modeller, ForceField, Simulation, PME, HBonds
from openmm.unit import nanometers, kelvin, picoseconds, femtoseconds
from openmm import Platform, LangevinIntegrator, Vec3, CustomExternalForce

Read PDB file:

pdb = PDBFile("./Chignolin.pdb")

Define the modeller from the PDB for adding H and solvent

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

forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
modeller.addHydrogens(pH=7.0)
modeller.addSolvent(forcefield=forcefield, neutralize=False, boxSize=Vec3(4.0, 4.0, 4.0)*nanometers)

PDBFile.writeFile(modeller.topology, modeller.positions, open("output.pdb", 'w'))

Create system

system = forcefield.createSystem(modeller.topology,
                                nonbondedMethod=PME,
                                nonbondedCutoff=1.0*nanometers,
                                constraints=HBonds
                               )

Check if forces uses periodic boundary conditions

system.getForce(2).usesPeriodicBoundaryConditions()

True

Define simulation parameters

integrator = LangevinIntegrator(300*kelvin, 1.0/picoseconds, 1.0*femtoseconds)

platform = Platform.getPlatformByName('CUDA')

simulation = Simulation(modeller.topology, system, integrator, platform)
simulation.context.setPositions(modeller.positions)

simulation.minimizeEnergy()

Run simulation

simulation.step(10_000_000)

Write outputs

positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(simulation.topology, positions, open("final_2.pdb", 'w'))

From the final_2.pdb file, it seems that all the water molecules have drifted way outside the 4 by 4 by 4 nm box initially used to solvation but I have no idea why this might be the case.

Hoping to get some help and advice on this...

chignolin.zip

[peastman]

See https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#periodic. It explains what's happening.