@gregoryross : Can you cut a release corresponding to the version of the code you used to generate results in the manuscript?

You'll also need to create a conda recipe with the versions of the tools (pinned) to what you used in the manuscript.

Most systems will require neutralizing counter ions, or may already have ions present. The topology of the ions in the Swapper class, which drives the MCMC exchanges, is the same as the topology as the water model. Currently, the ion parameters are hard-coded as the Joung and Cheatham models. To proceed, there are a couple of routes we could take depending on what we expect from users:

1. Users supply Swapper with a solvated system, possibly with salt already present. Swapper would then

• Read in the parameters of Na+ and Cl- from the supplied system (easy: this functionality is already possible, but not used) and
• Change the topologies of the salt ions to match the water model, and
• For each ion, add the positions of the hydrogens and any additional charge sites to the system coordinates.

2. Users supply Swapper with a solvated system without any Na+ and Cl- present. Swapper would then:

• Calculate the net charge of the system (easy)
• Convert the non-bonded parameters of enough water molecules to neutralize the system (also easy).

Given its simplicity, my preference is with option 2, but I’m open to suggestions. Option 2 has the benefit of ensures that the chemical potential is already calibrated for the ion parameters.

To make this repository easier to understand, navigate, and maintain, results and analysis (currently in the development folder) should be moved to a separate repository such that this repository only contains the source code, documentation, and unit tests.

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

I just realized this repo doesn't have travis testing.

@gregoryross : It would only take me ~20 min to get that set up for this repo. Should I take a stab at that now?

Can I get an update on the latest protocol recommendations (at least in the saltswap context) in terms of number of perturbation vs propagation kernels? As I understand it the issue is that perturbation remains slower than propagation, so there might be a sweet spot where there is still a lot of perturbation, but more propagation. We are exploring related issues for BLUES, and don't want to reinvent the wheel.

@sgill2 and I ran across a discussion of this in openmm/openmm#1832 (comment), where there was a graph of some (apparently older) data on how efficiency varied as a number of perturbation vs propagation kernels from @jchodera. Peter Eastman made the observation that it looked like what was important was the product of the two (though John noted that this was not quite the case). John mentioned that you guys would be exploring this more since perturbation is slower than propagation so there might be a sweet spot where one does more propagation than might otherwise be expected.

(cc @Limn1 )

We should add some basic sphinx docs on RTD following these guidelines from @bas-rustenburg.

Does this module still exist? Maybe it is for an older version? The examples don't seem to work without it.

@gregoryross : If the code is sufficiently stable for production calculations, we should at least cut an 0.1.0 release so we can conda-package this.

        atm_index = 0
        for atom in molecule:
>           charge = (1 - fraction) * initial_force[atm_index]["charge"] + fraction * final_force["charge"]
E           TypeError: list indices must be integers or slices, not str /home/rustenburg/miniconda3/lib/python3.6/site-packages/saltswap/swapper.py:811 

My guess is that it should be final_force[atm_index]["charge"] instead of final_force["charge"]. Will put in a PR with fix in a minute, and a unit test for this function to capture any other errors with this feature.

The README suggests the argument delta_chem to MCMCSampler has implicit units of kJ/mol. We should never have arguments carry implicit units---variables must always carry explicit units (using simtk.unit) or reduced units.

Trying the example run_sampler.py, i found two attributes that are not declared in the swapper class:
swapper.work_add
swapper.work_rm

run_sampler.py won't compile due to this issue.
This attributes are also present in NCMC Work Calculation.ipynb and SaltSwap Demo.ipynb.
Thanks

I{m having problems making the salt swapping work consistently. I have tried running the test systems provided (in testsystems), in particular DHFR at the same conditions shown in the published article. When running the script with a non-bonded cutoff distance of 10A and the following conditions: (c=100mM)
-c 0.1 -o out --ipdb minimized_dhfr.pdb --prmtop prmtop --water_name HOH -i 7500 -s 2000 -e 2000 --npert 1000 --nprop 10 --timestep 2.0 --save_freq 4 --platform CUDA
I get no swaping

If I change the cutoff distance to 10.1A, swapping is recovered but the distribution of concentrations and the correlation time are off when compared to the results published:

I'd like to mention that the system also fails to get swapping when using cutoff distances of: 12A,14A.
I'll attach all the files used.

dhfr.zip