Big jumps observed in extended particle coordinates, leading to NaN

[ljmartin]

Hi @craabreu ,
As discussed I've been running a few test systems. I came across a minimal example demonstrating a phenomenon where the extended particle coordinates occasionally jump a large distance (1 unit, being nanometers). Of course this leads to huge distortions in the system as it tries to catch up to the extended particle, which often causes a NaN.

The system is a charmm-gui lipid bilayer. The collective variables are simply x, y, and z coordinates of a sodium ion, made periodic since there is constant PBC. I track the value 1 integration step at a time, so at the bottom you can see that the extended particle coordinates jump first, shifting a large distance within the space of a single step.

The graph is in the notebok at https://github.com/ljmartin/tamd_jump , but here's a screen grab. Blue is the extended particle coordinates. The trajectory is also quite entertaining to watch on VMD or something.

The particular CV (x, y or z) that jumps seems to be stochastic, but in most cases you can see a huge jump in the extended particle, followed by the real CV trying to catch up.

And thanks for the library its been so useful !

[ljmartin]

I tried a couple of things to verify it's not something to do with my bad parameters:

• Reducing time step to 0.5 fs / step, to indicate whether or not the restraint is too tight leading to some kind of resonance
• Switching to a langevin integrator just in case
• Increasing the extended particle mass to 30,000 daltons to ensure it's not too light and getting "slingshotted" by the restraint
  But still see large jumps on the order of 1 nanometer.

[craabreu]

Hi @ljmartin

Thanks for reporting this issue. I wasn't able to run your example so far due to issues in my computer, but I spotted something that can lead to errors:

simulation.context.setPositions(state.getPositions())
simulation.context.setVelocitiesToTemperature(310*unit.kelvin, 1234)
simulation.context.setPeriodicBoxVectors(*state.getPeriodicBoxVectors())

First, there is no overriding for setPeriodicBoxVectors in _ExtendedSpaceContext, meaning that the global parameter Lx is not being changed as it should when you call this method. I just included this feature and pushed it to a different branch called fix_jumps.

Also, it is a good idea to call setPeriodicBoxVectors before calling setPositions because the latter uses the box vectors to set the values of the extended-particle coordinates. I will add a warning in the documentation.

BTW, this line should be corrected as well:
https://github.com/ljmartin/tamd_jump/blob/b9c32d2942a6fcfe71933254a2aecb289c4d1791/sodium_xyz.ipynb#L153

Would you please check if the error persists after all these changes?

Thanks.

[ajsilveira]

@ljmartin : I initially ran your script specifying first the box vectors and then the positions. I do this routinely, but I am not finding anything in the OpenMM documentation that says that box vectors should be specified before the positions. By doing this the NaN issue disappears, so this issue is independent of the jump observed in the extended particle coordinates. Regarding the jump in the coordinates, this issue is now fixed in the branch fix_jump, as @craabreu mentioned. This is the result I am getting:

[ljmartin]

Thanks @craabreu and @ajsilveira - I got the same result by:

• creating a System object from the (unequilibrated) PDB
• saving an equilibrated PDB, which writes the PBC in the first line
• recreate the system so UFEDMM can read the equilibrated box size from the pdb.topology during ufedmm.simulation instantiation
  Problem solved - thankyou!

While we're on this system, can I make a request for a setPositions(pos, extended=True) method? I noticed that the extended particle positions start a couple angstrom off the real CV. Extended particle=blue, real CV=orange:

Unless it's another issue with how I specified the box, I think the reason for the offset is the minimizer puts the extended particles into the minimum of an adjacent cell. I guess that should be equivalent, but notice that the zeroth extended particle should be starting at 4.05 but actually starts at 4.22 above.

PS - this only becomes an issue when I use 4fs timesteps, which I suppose is questionable with such tight restraints anyway!

[ajsilveira]

@ljmartin there are methods for specifying positions, velocities and biases (this is useful for restarting simulations). You just do:

ufed.set_positions(simulation, positions, extended=True)
ufed.set_random_velocities(simulation, seed)

data = np.load('bias.npy')
ufed.set_bias(simulation, data)

[ajsilveira]

@ljmartin Also, even though the simulation box is periodic, your collective variables are not. If your were computing the PMF for the translocation of the ion, for example, the positions at x=0 and x=L are not equivalent. In the context of ufed, metadynamics, etc., a periodic variable would be a torsion, for example, where 0 and 360 correspond exactly to the same state.

[craabreu]

@ljmartin and @ajsilveira, there is some confusion here due to changes in the API. The instructions for restarting positions, velocities, and biases won't work in the latest version. I'm sorry.

[craabreu]

Regarding the extended positions, they are not the same as the extended variable values. It is just a coincidence that, in this particular example, the collective variables are Cartesian coordinates as well. I will implement a method for recovering the extended variables from the context, so that they can be directly compared to the collective variables.

[craabreu]

In the meantime, a way of obtaining the extended-variable values from the extended positions is:

_, xcoords = simulation.context.getState(getPositions=True).getPositions(extended=True)
Lx = simulation.context.getParameter('Lx')
xvars = [v.evaluate(x, Lx) for v, x in zip(ufed.variables, xcoords)]

[craabreu]

Done! Now it is possible to do like this:

positions, xvars = simulation.context.getState(getPositions=True).getPositions(extended=True)
simulation.context.setPositions(positions, extended_positions=xvars)