Add alchemical ion support

[lohedges]

This PR adds support for adding alchemical ions to enforce charge neutrality during a perturbation. Currently the number of ions to apply to the perturbed state is configured via the --charge-difference command-line option. This specifies the difference in charge between the perturbed and reference states, i.e. charge(perturbed) minus charge(reference). For reference, the indices of the water molecules that are perturbed are written to the log. We also check the difference in end-state charge for all simulations and warn the user when there is a mismatch, i.e. suggesting to use --charge-difference if they want to preserve neutrality. In future, all of this could be automated if we want to completely remove the option, e.g. we could have some --charge-neutral flag that is True by default.

Currently we use the furthest N from X,Y,Z syntax to choose water molecules for perturbation. Here N is the number of ions that we need and X,Y,Z specifies the coordinates for the center of mass of the perturbable molecule. This results in a reproducible search for the same input. If we find that this causes issues, e.g. choosing a problematic water, then we could increase the number of molecules in the search, then pick randomly, e.g. first choosing 10N closest waters, then randomly choosing N from within this.

One point of difference with somd1 is the definition of the charge_difference option. Here it specifies the difference in charge between the perturbed and reference states, whereas for somd1 the difference appears to be reversed. (Which means the charge difference is the charge you need to add to the perturbed state to make it the same as the reference.) I'm happy to swap to whatever, but this definition seems less confusing to me and more in line with have the reference state actually be a reference.

The PR also contains a commit to save topology files for both end states to the output directory, i.e. system0.prm7 and system1.prm7. This allows a user to visualise or analyse trajectories using either topology as a reference. (This has been very useful for debugging angle distributions around perturbing atoms, for example.)

Happy for this PR to sit here so that @jmichel80 can comment before merging.

[lohedges]

Note that the test failures are expected:

• macOS is still using x86 for some reason, so there are no Sire or BioSimSpace packages available.
• Windows still fails on the test tear down due to the (unresolved) DCDFile unlinking issue. You can see that all tests pass before this occurs.

[mb2055]

From testing I can see a possible issue (unless I'm misunderstanding the functionality) - when using the --swap-end-states flag the stated charge difference appears to be the non-swapped difference.

For example, running the EG5 056~410 perturbation without the --charge-difference flag correctly raises the following warning: The charge difference of 1 between the end states does not match the expected value of 0. Please specify the 'charge_difference' if you wish to maintain charge neutrality. .
Adding the --swap-end-states flag and running again without the --charge-difference flag raises the exact same warning The charge difference of 1 between the end states does not match the expected value of 0. Please specify the 'charge_difference' if you wish to maintain charge neutrality. Should this not be The charge difference of -1...? Or is the alchemical ion also included in the swapped end states?

[lohedges]

I think this is okay. The swapping is only done internally by Sire dynamics itself, i.e. the input system is not swapped. The flag is simply used to swap the direction of the perturbation that is passed in. We should consider all of the setup stage, i.e. pre dynamics to operate on the system in the input order, i.e. non-swapped.

[mb2055]

Okay, that makes sense.
So, just for the sake of clarity, the --swap-end-states flag applies to all perturbable molecules in the system, including any co-alchemical ions?

[lohedges]

Yes, exactly. The Sire code doesn't care about how many perturbable molecules there are. When setting them up in converting to a OpenMMMolecule it will just invert the properties (other than coordintes) if swap_end_states is True.

[mb2055]

Good to go 👍