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')>})})
I've updated the example to reflect changes in PR openforcefield/openff-interchange#30
PotentialHandler
object now operate on a correspondingParameterHandler
instead of the entireForceField
ForceField.create_openff_system
method to demonstrate one potential home for this functionality