Fix issue #86

[lohedges]

This PR closes #86.

This is a first pass to handle molecular coordinates after vacuum simulation with GROMACS. We implement vacuum simulations by using pseudo-PBC conditions, i.e. running calculations in a near infinite box, as described here, which is works on all version of GROMACS that we support. However, absolute molecular coordinates can end up being very large after simulation, which causes problems when trying to write them to certain formats, i.e. the fixed-width formatting restrictions are violated.

In this PR we compute the center of mass of the system following the simulation, then translate this back to the center of the box. I've needed to take care in the Exscientia Sandpit since they use two approaches read molecular frames for getSystem. One uses the trajectory file with wrapped coordinates, the other uses the final gro frame, in which the coordinates are unwrapped. As such, one center of mass calculation needs to use a periodic space, the other does not.

Issues that I anticipate with this approach that we might need to revisit at a later date are:

• For multi-molecule systems a center of mass translation still might end up with large absolute coordinates, e.g. if the molecules aren't close to each other. I doubt this is the typical use case for vacuum simulations, so we'll see if it becomes an issue. Here we might be able to translate molecules independently, forcing some minimum distance between them afterwards.
• This fix doesn't apply for trajectory frames obtained via getTrajectory. Something similar could be done here too, but the user would lose the ability to track displacements along a trajectory, etc.
• This might not work if the user sets some config option to change the wrapping during simulation. But with customised simulations a lot of our functionality is a bit hard to guarantee.

Ultimately we might want to expose a centerMolecules function or similar, which the user can apply before trying to write to file. This could even be an option within saveMolecules? Let me know if you think this would be a better solution, i.e. the user would need to opt-in, rather than us modify something. In this case I feel like they should just expect things to work, though, and it would be annoying if the behaviour was different to other engines. (I'm not sure how they deal with molecules wandering off during vacuum simulations. Perhaps they internally recenter based on the center of mass periodically?)

The center of mass calculation is currently written in Python, but can take advantage of the C++ code once it is implemented there. I've also updated the funnel metadynamics setup code to use the new functionality so as to avoid code duplication. This is tested via tutorial notebook, which I've run independently.

Checklist

• I confirm that I have merged the latest version of devel into this branch before issuing this pull request (e.g. by running git pull origin devel): [y]
• I confirm that I have permission to release this code under the GPL3 license: [y]

Suggested reviewers:

@chryswoods

[chryswoods]

Thanks - I agree this should eventually be more general, as I've seen this myself in sire/openmm gas phase simulations I've run. I agree that something that periodically translates back to center is the least worst solution.

[lohedges]

Good to know that this isn't isolated to GROMACS. I'll think of a more general solution if it crops up again. Bizarrely, for the test molecule here, the translation to the box edge happens within 50 steps of minimisation. It's not blown up, so I'm surprised that it's moved so far so quickly when no dynamics are being performed.

[chryswoods]

That's really weird. Minimisation shouldn't give the molecule any linear momentum as it isn't sat in an external field? Certainly not enough to pull it so far outside the box. Do we see the same thing with mols.minimisation().run()? If not, then this may be worth raising with the gromacs developers as a bug.

[lohedges]

I think the issue here was that the default box center for GROMACS is (x/2, y/2, z/2), so I should have probably translated the system to the center before the simulation, then back again afterwards. That said, the final coordinates that I was seeing would have still implied that it had moved a long way, rather than just drifting off the edge of the periodic box. I'll test both approaches before deciding what is the best way of doing things.