Unstable saltswap simulations in octahedral solvent boxes

[matiasmachado]

Hi! I've been testing saltswap in systems from the reference paper and some of my own interest.

In doing so, I've found that using octahedral solvent boxes instead of rectangular ones introduces some kind of instability during the simulations.

No matter the amount of equilibration steps nor the softness of the protocol, the perturbation or propagation steps, all systems sooner or later end up crashing with the following error message:

Traceback (most recent call last):

  File ".conda/envs/openmm/lib/python3.6/site-packages/saltswap-0.5.2.1-py3.6.egg/saltswap/swapper.py", line 656, in attempt_identity_swap  work_measurement=self.work_measurement)
  File ".conda/envs/openmm/lib/python3.6/site-packages/saltswap-0.5.2.1-py3.6.egg/saltswap/swapper.py", line 536, in _ncmc
    self.ncmc_integrator.step(nprop)
  File ".conda/envs/openmm/lib/python3.6/site-packages/simtk/openmm/openmm.py", line 15792, in step
    return _openmm.CustomIntegrator_step(self, steps)
Exception: Particle coordinate is nan

During handling of the above exception, another exception occurred:
Traceback (most recent call last):
  File "sample_amber_system_v2.py", line 187, in <module>
    salinator.update(nattempts=1)
  File ".conda/envs/openmm/lib/python3.6/site-packages/saltswap-0.5.2.1-py3.6.egg/saltswap/wrappers.py", line 405, in update
    self.swapper.update(self.context, nattempts=nattempts, cost=chemical_potential, saltmax=saltmax)
  File ".conda/envs/openmm/lib/python3.6/site-packages/saltswap-0.5.2.1-py3.6.egg/saltswap/swapper.py", line 915, in update
    self.attempt_identity_swap(context, penalty=cost, saltmax=saltmax)
  File ".conda/envs/openmm/lib/python3.6/site-packages/saltswap-0.5.2.1-py3.6.egg/saltswap/swapper.py", line 661, in attempt_identity_swap
    if detail[0] == 'Particle coordinate is nan':
TypeError: 'Exception' object is not subscriptable

[1]+  Exit 1                  python sample_amber_system_v2.py -c 0.2 -o out --prmtop DHFR_octbox.prmtop --inpcrd DHFR_octbox.inpcrd --water_name HOH -i 7500 -s 2000 -e 2000 --npert 1000 --nprop 10 --timestep 2.0 --save_freq 4 --platform CUDA

Despite the number of iterations to reach the error depends on each particular run, it trends to be lower than 1000 (i.e. 4ns).

I've also tried different GPU platforms ( GTX 1080, Titan Xp) but got the same results.

I've created the following crashing test corresponding to the DHFR system in a octahedral solvent box with almost the same water amount of the "minimized_dhfr.pdb". I've just modified sample_amber_system.py to record the iteration steps in a .log file.
octbox.zip

I'm using openmm 7.4.0 with openmmtools 0.19.0 (see my environment.yml for more details)
environment.zip

I also patched saltswap/swapper.py (version 0.5.2) with self.ncmc_integrator.reset() according to #29