Skip to content

Instantly share code, notes, and snippets.

@mattwthompson
Last active October 21, 2020 22:55
Show Gist options
  • Save mattwthompson/c603d4895d9eafcc37fe7ab8afcb64d8 to your computer and use it in GitHub Desktop.
Save mattwthompson/c603d4895d9eafcc37fe7ab8afcb64d8 to your computer and use it in GitHub Desktop.
  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')>})})
@mattwthompson
Copy link
Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment