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')>})})
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
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)]
andmy_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.