Adding a new VirtualSite definition

[Tristan Croll]

I'd like to look into adding a new type of VirtualSite to provide a strict symmetry implementation. My (perhaps not fully-formed) idea is that each virtual particle would have a position defined by the symmetry-transformed coordinates of its parent atom, and act to collect all the extra nonbonded interactions with its symmetry neighbours. The concept is fairly straightforward (barring a few wrinkles in dealing with locations where the symmetry breaks), but I've looked at enough of the code to see that the implementation would require edits in quite a few different places. Just wondering if someone more familiar with the architecture could give me a few pointers to get started?

[Peter Eastman]

That sounds interesting! Can you give more details on how it would work? What sorts of symmetry do you have in mind?

To implement a new type of virtual site, there are two main elements you need to provide. The first one is the formula for computing the virtual site's position based on the positions of a set of other atoms. So let's start with that. How will the virtual site position be calculated?

The second element, which is often more complicated, is that you need to invert the above formula so you can take forces acting on the virtual site, and properly distribute them back to the atoms it was computed from. In principle that's straightforward once you've decided on the formula, but in practice it can end up being pretty complicated.

[Tristan Croll]

The purpose is to cleanly handle symmetry interfaces in crystal structures (and EM structures with strict non-crystallographic symmetry). My particular focus is on using interactive MD for model building/refinement into experimental maps, but I'm sure there would be other uses for it. In the model-building world, you want to be working with a single asymmetric unit (the minimum repeat unit necessary to replicate the entire structure after applying symmetry operators), while still correctly handling the sites where it contacts symmetry neighbours. I thought about trying to do this with strong cartesian symmetry restraints, but the more I think about it the more I think this approach will work out much better.

The position formula and its inversion is fairly straightforward - just a rigid-body transform of the parent atom coordinate (3x3 rotation matrix and 1x3 translation) and its inverse. All virtual atoms in one symmetry group would of course have the same transform, and would not interact with each other - just with "real" atoms and virtual atoms from other symmetry groups. I think this can be handled using CustomNonbondedForce interaction groups?

[Tristan Croll]

This video might help: https://drive.google.com/open?id=1d9QLl ... ztVRNoodjt. It's a screen capture of what I have working right now. It's a live running simulation, where the bright coloured atoms are the "real" molecule and the darker ones are the symmetry equivalents. At the moment the symmetry atoms are shown for information only, but to make this as good as it can be I need to move beyond that.

[Peter Eastman]

I see what you mean. We actually discussed something along these lines a while back, but never got as far as implementing it. See https://github.com/pandegroup/openmm/issues/1437.

There are a lot of pieces that will be needed eventually, but here's where to start.

1. Modify VirtualSite.h and VirtualSite.cpp to add the new class to represent the new kind of virtual site.

2. Add the reference implementation of it to ReferenceVirtualSites.cpp. You'll add a new block in computePositions(), and another one in distributeForces().

3. Update TestVirtualSites.h to test your code and make sure it works. It includes one method for every virtual site type that creates a simple system with just a few particles and makes sure all the results are correct. You'll add another one for the new virtual site type. You should also update testConservationLaws(). It creates a system with one of every virtual site type, then runs a simulation and checks that energy, momentum, and angular moment are all conserved. Bugs in distributing forces will usually cause that to fail. So add another molecule with your new virtual site type.

That should be plenty to start. Once that's working, you can move on to the other platforms, serialization, etc. (See the checklist at https://github.com/pandegroup/openmm/wi ... ew-feature for everything that goes into a new feature.)

[Tristan Croll]

Ha - that's exactly what I was talking about!

To add a little to the discussion from your GitHub link, there are all sorts of extra complications involved in taking a deposited structure and reconstructing a full P1 box in this way. In no particular order:

- it's not all that uncommon for a deposited structure to have one or more residues with broken symmetry at the interface with symmetry copies. A classic example is where an asparagine or glutamine side chain H-bonds to its own rotationally symmetric equivalent. In that case one terminal amide has to be flipped 180 degrees relative to the other. Interestingly enough, the PDB architecture doesn't have the framework to annotate such cases, so you can't even rely on reading the metadata to find them. There are various other similar situations, where you'd either need to start implementing multiple non-interacting conformers of the offending residues, or doing a preliminary expansion to a lower space group, assignment of compatible rotamers, and then work up. Somewhat painful either way.

- It's also (somewhat more) common to find symmetric small molecules/ions on special positions, where one or more symmetry operators maps the position back onto itself. Again, something that needs to be caught and specially handled.

Nothing that's impossible to handle, but lots of tricky subtleties.

Anyway, thanks for the pointers. Will start looking into it in earnest hopefully in a few weeks time.

[Peter Eastman]

Thanks for the helpful insights! If you feel inclined to post those comments on the github issue as well, that would be really helpful. That way everything will be in one place.