Jumping Molecule

[Alexander Falk]

Hello All,
I ran a simulation of long protein chain in a water box. On inspecting the simulation results, I noticed that the chain seems to jump between two frames (picture attached), at about the 30 ns mark of a 225 ns simulation. The simulation seemed to run smoothly and the output looked normal otherwise, but I am trying to figure out what happened here and undo this jump (if possible) because I think it will get in the way of my analysis.
I know these sorts of jumps can happen due to the PBC, but I do not know if that is what is going on here, because the jump size is much smaller than the PBC box and the protein appears to have been rotated as well. Also, I tried to resolve this using the trajMD and the VMD pbctools plugin, but was not able to do so with either.
At the same time, I do not think this is actually just the protein moving in the simulation, the time difference between these frames is only 20 ps and the simulation is at normal temperature and pressure. Further, there are positional restraints on some of the atoms that have been moved, and they never should have been allowed to move that far.
If anyone can help me figure out what happened/how to fix it, I would greatly appreciate it. Sorry if I was unclear about anything. I can provide more details if needed.
Thank you,
Sandy

[Peter Eastman]

I'm kind of mystified by this. As you say, periodic boundary conditions should only produce translations, not rotations. Also, you said you included position restraints. How are those restraints implemented? What's the exact form of them? Can you evaluate the energy of the restraints before and after it moved?

[Alexander Falk]

Hi Peter,
Thank you so much for the follow up. So I actually just figured out what seems to have been the issue (although I had to rerun the sim). I solvated the system using the command:

modeller.addSolvent(forcefield, boxSize=Vec3(7.0, 7.0, 21.0)*nanometers, model='tip4pew',neutralize=True)

My protein started the simulation sitting right along the PBC edge (essentially stretched from (0,0,10) to (0,0,100) in an extended conformation). When I noticed this puts the protein right on the edges of the box, I guessed it may be confusing the system (although it should not be an issue in theory?). So I adjusted the protein starting position so that it was stretched from (35,35,10) to (35,35,100) and reran the sim, with otherwise the exact same parameters. This seems to have worked and I do not see any "jumping" in the trajectory (the protein otherwise acts the same as it did in the previous sim).

To address your questions though, I implemented my restraints using the following code:

force = openmm.CustomExternalForce("k*((x-x0)^2+(y-y0)^2+(z-z0)^2)")
force.addGlobalParameter("k", .5*kilocalories_per_mole/angstroms**2)
force.addPerParticleParameter("x0")
force.addPerParticleParameter("y0")
force.addPerParticleParameter("z0")
for i, atom_crd in enumerate(modeller.positions):
if (i == 1 or i == 5 or i == 18):
force.addParticle(i, atom_crd.value_in_unit(nanometers))
system.addForce(force)

I am not sure if this is the proper way to do it, but outside of the "jumping" these restraints seem to be effective. I am not sure how to pull energies for the restraints in particular (is there a reporter for this?), but I did not see any energy change in the overall simulation corresponding to the jump. Also, restraints aside, the amount the protein moves in the pre and post-jump frames as way more than the movement seem between any other two adjacent frames.

[Peter Eastman]

Glad you got it working! That restraint force is going to cause problems with periodic boundary conditions. Instead define it as

Code: Select all

force = CustomExternalForce("k*periodicdistance(x, y, z, x0, y0, z0)^2")

There's a discussion of this at http://docs.openmm.org/latest/api-pytho ... ernalForce.

[Alexander Falk]

Awesome, thank you.