Add link-atom method for molecules crossing ML-MM boundary

[epretti]

Initial implementation of the link-atom method. Not ready; documentation to be added once API is finalized.

[peastman]

The way extra atoms are added to the Topology is going to cause problems. The atoms really need to be sorted by residue. You keep the list of existing atoms as they are, then add more atoms at the end of the list belonging to many different residues. Think of what would happen if you saved that Topology to a PDB file, for example. It would come out in the wrong order and wouldn't be interpreted correctly.

We should allow returning the new Topology, for example by passing returnTopology=True. People will also need to update atom positions to be correctly ordered for the new Topology. Perhaps return a mapping from old atom indices to new atom indices?

[epretti]

For returning a new Topology and index mapping, perhaps something like returnInfo=True to return a {'system': ..., 'topology': ..., 'oldToNew': ...} instead of just a System? That could make it less painful for other embeddings to return different kinds of extra information in the future.

I've now realized that inserting virtual sites into each residue in order isn't straightforward. We have to add the virtual sites to the System too, which is easy when they are appended to the end but difficult otherwise, since the System is user-provided. A messy and limited approach might be to process each Force, remapping all of the atom indices, and giving up if we see an unknown kind of Force. Alternatively, we could add a new "chain" to the end of the Topology with just the virtual sites, which is technically far easier and more flexible, although it adds meaningless items to the Topology. I would probably prefer the latter, unless there is an easier way to accomplish the former, or a better approach in general I haven't thought of yet.

[epretti]

There have been a few subtleties that I have run into while testing this.

First, configurations sampled by a system using the link-atom method may differ from those sampled from a standard simulation even if the standard simulation, link atom ML subset, and link atom MM subset all use the same force field, since the link atoms may interact with the rest of the ML subset: even though the distance between a link atom and its ML-side link bond atom is fixed, distances to the other ML subset atoms are not. This seems intrinsic to the method and shouldn't normally be a major problem; it just took me a bit of time to ensure that there was not a bug in the implementation.

Second, checking for which bonds to remove from the MM force field based on the ML subset + MM-side link bond atoms may be incorrect in some cases.

Consider a torsion like $-A_{ML}-B_{MM}-C_{MM}-D_{ML}-$. The current rule would remove this torsion from the ML subset even though it wouldn't be evaluated by the ML potential (which would see separate pieces $-A_{ML}-H$ and $H-D_{ML}-$). The rule needs to be changed to correctly account for which interactions are actually part of the ML region.

The angle $A_{ML}-B_{MM}-C_{ML}$ suffers from the same problem. However, I think it may cause a more severe issue, since the link atoms will be constrained to be close to each other, which could distort the molecule. I wonder if we should warn or raise an error for such cases (multiple link bonds with the same MM-side atom)?

[peastman]

I think the first problem is intrinsic to the method. The ML potential is seeing a capped molecule, which is not the same as the original. We hope that it's close, but it's not the same.

For the second problem, I think the rule is that an angle or torsion should only be removed if it will be calculated by the ML potential. That means it consists exclusively of ML atoms. A capping hydrogen can take the place of an MM atom, but only at the end of the sequence. In your example $-A_{ML}-B_{MM}-C_{MM}-D_{ML}-$, B and C would be replaced by hydrogens in the ML region, but there would be no bond between them. Therefore the ML potential would not calculate that torsion. But if you change it to $-A_{MM}-B_{ML}-C_{ML}-D_{MM}-$, that would get computed by ML. A and D would be replaced by hydrogens and the full torsion A-B-C-D would be present in the ML region. So the MM atoms have to be at the ends, not in the middle.

However, I think it may cause a more severe issue, since the link atoms will be constrained to be close to each other, which could distort the molecule.

Yes, I think it's reasonable to require that a MM atom can only be bonded to a single ML atom.

[epretti]

I think the rule is that an angle or torsion should only be removed if it will be calculated by the ML potential.

Agreed. I thought so too but failed to consider the cases where my first attempt does something different. This turns out to be tricky since one needs to identify improper torsions also. The implementation now checks explicitly for these; I tried to think if there was a simpler correct way to do this but I don't know if there is.

This should be ready for further review.

[peastman]

I see what you mean about impropers, especially that they have to be identified from the Topology, not the System.

[peastman]

We should handle CustomNonbondedForce too, since MechanicalEmbedding supports it. This likely requires adding exclusions to prevent the link atom from interacting with anything else.

[peastman]

This implies that some implementations might modify the Topology. We should make it clear that they should only create a new Topology, not modify the existing one.

[epretti]

I've made changes to support CustomNonbondedForce. I also added a test for it in the ordinary (non-link-atom) case since there wasn't any. It turned out that there was a typo in a method call for this case, which should now be fixed.

[peastman]

Looks good. I wonder if we should always use the exclusions, even if it's only a NonbondedForce? Setting the parameters to 0 will prevent it from interacting with anything, unless the distance to another particle becomes exactly 0, in which case the energy becomes nan. Perhaps that wouldn't actually happen in practice thanks to the repulsion from the parent atom?

[epretti]

I would imagine that shouldn't be possible because the site should be sandwiched between two real atoms. I considered doing this but was worried about performance. It seems like adding an exception from one particle to every other can slow down OpenMM significantly. (For instance, adding an exception from one particle to all others in a 12,255 atom water box takes OpenMM from about 1270 to 530 ns/day. I am curious as to why and if that can be mitigated, but that's an entirely different issue.)