minimal_handler_example.md

1. Import OpenFF objects
2. Populate a SMIRNOFFBondHandler
3. Store the output of running typing
4. Store the actual potential objects

In [1]: from openforcefield.typing.engines.smirnoff import ForceField
^[[A^[[A^[[A^[[A^[[A

In [2]: from openforcefield.topology import Topology, Molecule

In [3]: parsley = ForceField('openff-1.2.1.offxml')

In [4]: top = Topology.from_molecules([Molecule.from_smiles('CCO'), Molecule.from_smiles('CC')])

In [5]: from openff.system.components.smirnoff import SMIRNOFFBondHandler

In [6]: bonds = SMIRNOFFBondHandler()  # This class inherits from a common PotentialHandler class

In [7]: bonds.store_matches(parameter_handler=parsley['Bonds'], topology=top)  # All subclasses must implement this method

In [8]: bonds.slot_map  # key-val pairs of unique identifiers for bonds ("slots") and potentials (SMIRKS, in SMIRNOFF world) 
Out[8]:
{(0, 1): '[#6X4:1]-[#6X4:2]',
 (0, 3): '[#6X4:1]-[#1:2]',
 (0, 4): '[#6X4:1]-[#1:2]',
 (0, 5): '[#6X4:1]-[#1:2]',
 (1, 2): '[#6:1]-[#8:2]',
 (1, 6): '[#6X4:1]-[#1:2]',
 (1, 7): '[#6X4:1]-[#1:2]',
 (2, 8): '[#8:1]-[#1:2]',
 (9, 10): '[#6X4:1]-[#6X4:2]',
 (9, 11): '[#6X4:1]-[#1:2]',
 (9, 12): '[#6X4:1]-[#1:2]',
 (9, 13): '[#6X4:1]-[#1:2]',
 (10, 14): '[#6X4:1]-[#1:2]',
 (10, 15): '[#6X4:1]-[#1:2]',
 (10, 16): '[#6X4:1]-[#1:2]'}

In [9]: bonds.store_potentials(parameter_handler=parsley['Bonds'])  # and also must implement this method ...

In [10]: bonds.potentials  # ... which serves as a mapping between the potential identifier and the parametrized force field data
Out[10]:
{'[#6X4:1]-[#6X4:2]': Potential(parameters={'k': <Quantity(526.418876, 'kilocalorie / angstrom ** 2 / mole')>, 'length': <Quantity(1.52108706, 'angstrom')>}),
 '[#6X4:1]-[#1:2]': Potential(parameters={'k': <Quantity(753.788152, 'kilocalorie / angstrom ** 2 / mole')>, 'length': <Quantity(1.09422343, 'angstrom')>}),
 '[#6:1]-[#8:2]': Potential(parameters={'k': <Quantity(660.456413, 'kilocalorie / angstrom ** 2 / mole')>, 'length': <Quantity(1.42841554, 'angstrom')>}),
 '[#8:1]-[#1:2]': Potential(parameters={'k': <Quantity(1113.99071, 'kilocalorie / angstrom ** 2 / mole')>, 'length': <Quantity(0.973822581, 'angstrom')>})}

In [11]: bonds.potentials[bonds.slot_map[(2, 8)]].parameters['k']  # Look up the force constant on some known bond (which itself may not know its physics!)
Out[11]: 1113.990711326 <Unit('kilocalorie / angstrom ** 2 / mole')>

In [12]: from openff.system.stubs import *  # Monkey-patch new methods onto existing OpenFF objects

In [13]: system = parsley.create_openff_system(topology=top)

In [14]: system.handlers['Bonds']  # everything that this PotentialHandler knows
Out[14]: SMIRNOFFBondHandler(name='Bonds', expression='1/2 * k * (r - length) ** 2', independent_variables={'r'}, slot_map={(0, 1): '[#6X4:1]-[#6X4:2]', (0, 3): '[#6X4:1]-[#1:2]', (0, 4): '[#6X4:1]-[#1:2]', (0, 5): '[#6X4:1]-[#1:2]', (1, 2): '[#6:1]-[#8:2]', (1, 6): '[#6X4:1]-[#1:2]', (1, 7): '[#6X4:1]-[#1:2]', (2, 8): '[#8:1]-[#1:2]', (9, 10): '[#6X4:1]-[#6X4:2]', (9, 11): '[#6X4:1]-[#1:2]', (9, 12): '[#6X4:1]-[#1:2]', (9, 13): '[#6X4:1]-[#1:2]', (10, 14): '[#6X4:1]-[#1:2]', (10, 15): '[#6X4:1]-[#1:2]', (10, 16): '[#6X4:1]-[#1:2]'}, potentials={'[#6X4:1]-[#6X4:2]': Potential(parameters={'k': <Quantity(526.418876, 'kilocalorie / angstrom ** 2 / mole')>, 'length': <Quantity(1.52108706, 'angstrom')>}), '[#6X4:1]-[#1:2]': Potential(parameters={'k': <Quantity(753.788152, 'kilocalorie / angstrom ** 2 / mole')>, 'length': <Quantity(1.09422343, 'angstrom')>}), '[#6:1]-[#8:2]': Potential(parameters={'k': <Quantity(660.456413, 'kilocalorie / angstrom ** 2 / mole')>, 'length': <Quantity(1.42841554, 'angstrom')>}), '[#8:1]-[#1:2]': Potential(parameters={'k': <Quantity(1113.99071, 'kilocalorie / angstrom ** 2 / mole')>, 'length': <Quantity(0.973822581, 'angstrom')>})})

Regarding those "slots" objects (the keys like (0, 1), would it be helpful for them to encode something about "what to do" with those particle indices? This would help distinguish bonds/"pairwise distance" slots from, say, a nonbonded exception, LibraryCharge, or ChargeIncrement that also applied to (0, 1). Distinguishing this might make it easier to have meaningful "key collisions", and they could potentially also encode the geometric transform between coordinates and the measurement of interest (this could help distinguish between dihedral-angle and out-of-plane style impropers)

My current opinion is that I'd prefer to keep the topology objects blissfully unaware of their physics until there is a compelling reason to add in that complexity

Avoiding key collisions

My hope - that has not really been stress-tested through implementation - is that there will be no high-level tracking of slots, instead only tracking them within a handler. This example is already restricted to the context of handler, so it's not really demonstrated here. In this case, a bond slot (tuple of atom indices) can only be accessed through the handler (i.e. my_system.handlers['Bonds'].slot_map[(2, 8)]. So, even when keys in one handler are similar or identical to other keys in other handlers, their intended use is implied by the handler they're inside of. Similar or identical "slot keys" can exist, but only in separated contexts, i.e. my_system.handlers['Bonds'].slot_map[(2, 8)] and my_system.handlers['ChargeIncrement'].slot_map[(2, 8)] will not know about each other, and if a high-level lookup is added, it will need to know both a handler and slot.

Getting away with this may be wishful thinking - happy to discuss breaking/problematic counter-examples here.

I've updated the example to reflect changes in PR openforcefield/openff-interchange#30

• The topology is now two molecules (an ethane and an ethanol)
• The methods of the PotentialHandler object now operate on a corresponding ParameterHandler instead of the entire ForceField
• Monkey-patched a ForceField.create_openff_system method to demonstrate one potential home for this functionality