VirtualSite support

[trevorgokey]

Checklists

• Add tests

• Update docstrings/documentation, if applicable

• Update changelog

• Update the SMIRNOFF specification

• Implement/parse the VirtualSite offxml specification

• Implement create_force

• Fixes Orphaned TopologyParticle? #651

This PR is very experimental and can change without warning.
Edit: This PR is near completion, so functionality will likely not be changed when merged.

• Ready to go!

Old/Outdated notes:

Limitations noticed:

It seems to not be possible to apply two virtual sites using the same virtual site type and smirks combination. For example, my cursory testing shows that, for a triplet of atoms (e.g. a water molecule), it is only possible to apply exactly one DivalentLonePair virtual site. This limits the ability to add multiple virtual sites on a single pattern.

Internal representation:

The current implementation is represented internally, and to the user, as the following. This is applying random VirtualSites to a water molecule, and below is what what is returned by label_molecules for VirtualSites only:

*** VirtualSites
{'store': OrderedDict([((0, 1),
                        <VirtualSiteBondChargeType with smirks: [O:1]-[H:2]  epsilon: 0.05 kcal/mol  distance: 0.3 A  chargeincrement1: 0.2 e  chargeincrement2: 0.1 e  type: BondCharge  sigma: 0.1 A  >),
                       ((0, 1, 2),
                        <VirtualSiteDivalentLonePairType with smirks: [H:1]-[O:2]-[H:3]  epsilon: 0.05 kcal/mol  distance: 0.3 A  chargeincrement1: 0.2 e  chargeincrement2: 0.1 e  type: DivalentLonePair  outOfPlaneAngle: 1.0 deg  inPlaneAngle: 0.0 deg  sigma: 0.1 A  >),
                       ((1, 2),
                        <VirtualSiteBondChargeType with smirks: [O:1]-[H:2]  epsilon: 0.05 kcal/mol  distance: 0.3 A  chargeincrement1: 0.2 e  chargeincrement2: 0.1 e  type: BondCharge  sigma: 0.1 A  >)])}

The spec. shows elements as simply VirtualSite, but I discovered that I needed to specialize the elements based on the type relatively early for validation purposes. This differs somewhat what the specification/xml represents, so I am not sure if the above representation is desirable.

[lgtm-com]

This pull request introduces 1 alert when merging 96f76a6 into 20712c1 - view on LGTM.com

new alerts:

• 1 for Unused local variable

[codecov-io]

Codecov Report

Merging #548 into master will not change coverage by %.
The diff coverage is n/a.

[j-wags]

Awesome progress. Recapping our meeting today:

• We should allow multiple matches for VirtualSites, and multiple of the exact same SMIRKS
• When initializing a ForceField, should we deduplicate VirtualSites with same SMIRKS, type, angles, distances, chargeincrements?
  • We should not deduplicate VirtualSites during VirtualSiteType or VirtualSiteHandler initialization
• Does OpenMM have a problem if two virtualsites are really close to each other (eg. do they cause an energy explosion?) Under what conditions does OpenMM automatically deduplicate or apply nonbonded exceptions between virtualsites?
  • Unresolved -- We should check
• Should we throw an exception if two virtualsites are defined on top of each other (same reference atoms, angles, distances) during system creation?
  • Yes, but add kwarg to create_openmm_system for allow_overlapping_virtual_sites
• We should change the tagname of the “VirtualSite” elements to be distinct -- that is, the XML element should be titled “BondTypeVirtualSite” etc. We should change the SMIRNOFF spec docs accordingly.

I had some thoughts on implementation after our meeting:

• If at all possible, the contents of the VirtualSiteType subclasses should contain anything essential about the functional form of the VirtualSite. That way, if a user wants to experiment by making a new VirtualSite subclass, they shouldn’t need to change OpenFF toolkit code, but rather they should be able to “register” a new VirtualSite subclass, and then the Toolkit should automatically load it when initializing a ForceField and recognize the new unique tagname.
  • I don’t know exactly what the OpenMM API for defining VirtualSites looks like, but it would be ideal if each VirtualSiteType subclass contained the necessary OpenMM keywords/classes/methods to add the proper VirtualSite to an OpenMM System. Right now, we kinda do this with the ParameterHandler _OPENMMTYPE class attribute, where we just hold on to a “type” object for each OpenMM Force, and repeatedly throw keywords at its constructor during system creation for each SMIRKS match we find.
  • We may need a more sophisticated method for doing this depending on how the OpenMM VirtualSite API looks. If it’s really hairy, we might want VirtualSiteType to be the first ParameterTYPE with its own create_force method

Very roughly, an API example could look like this:

from openforcefield.topology import Molecule
from openforcefield.typing.engines.smirnoff import ForceField 

class MyNewVirtualSiteHandler(VirtualSiteHandler):
    _TAGNAME = ”VirtualSites”
    class MyNewVirtualSiteType(VirtualSiteType):
        _ELEMENT_NAME = ”MyNewVirtualSite”
        def _create_force(self, openmm_system, particle_indices):
            (Here the OpenMM System is modified to contain the desired virtual site)

my_topology = Molecule.from_smiles(‘HOH’).to_topology()
my_ff = ForceField(‘my_ff_with_my_new_virtualsite_type.offxml’, parameter_handler_classes=[BondHandler, AngleHandler… MyNewVirtualSiteHandler])
system = my_ff.create_openmm_system(my_topology)

• To enable a plugin-like architecture for new VirtualSite types, I think that populating a dictionary/registry with available VirtualSite subclasses would work. We do something like this when initializing ParameterHandlers. This mapping could have keys of VirtualSite tagnames and values of VirtualSiteType subclasses. Here are some existing code snippets that do something similar:
  • “Registering” all parameterhandlers
  • Lookup for parameterhandlers
  • Finding all subclasses (linking just for reference, we could just call this function)

[j-wags]

LGTM, @trevorgokey

...

On this pale fear seized every one; they were so frightened that their arms dropped from their hands and fell upon the ground at the sound of the goddess's voice, and they fled back to the city for their lives. But Ulysses gave a great cry, and gathering himself together swooped down like a soaring eagle. Then the son of Saturn sent a thunderbolt of fire that fell just in front of Minerva, so she said to Ulysses, "Ulysses, noble son of Laertes, stop this warful strife, or Jove will be angry with you."

Thus spoke Minerva, and Ulysses obeyed her gladly. Then Minerva assumed the form and voice of Mentor, and presently made a covenant of peace between the two contending parties.

THE END

[trevorgokey]

Merged! Amazing. Thanks @j-wags and @mattwthompson for all the discussions and feedback!