The Open Force Field toolkit applies parameters based on the SMARTS/SMIRKS language. These examples are intended to expose how ForceField assigns SMIRKS-based parameters to test molecules.
This example shows the parameters applied to a single molecule that is created from a SMILES string.
The first output section consists of Bond parameters that would be assigned to the molecule. Each bond is described by the indices of the two atoms it connects, its "parameter ID", and the SMIRKS that caused it to match the target bond.
The other ParameterHandlers print in order after that.
Note that the ToolkitAM1BCC and Electrostatics handlers do not assign SMIRKS-based parameters, therefore they do not print any match information.
In [1]:
from openforcefield.topology import Molecule, Topology
from openforcefield.typing.engines.smirnoff import ForceField
# Create a simple molecule from SMILES and turn it into a topology.
molecule = Molecule.from_smiles('C=C(N)(C)')
topology = Topology.from_molecules([molecule])
# Let's label using the smirnoff99Frosst force field
forcefield = ForceField('openff-1.0.0.offxml')
# Run the molecule labeling
molecule_force_list = forcefield.label_molecules(topology)
# Print out a formatted description of the parameters applied to this molecule
for mol_idx, mol_forces in enumerate(molecule_force_list):
print(f'Forces for molecule {mol_idx}')
for force_tag, force_dict in mol_forces.items():
print(f"\n{force_tag}:")
for (atom_indices, parameter) in force_dict.items():
atomstr=''
for idx in atom_indices:
atomstr += '%3s' % idx
print("atoms: %s parameter_id: %s smirks %s" % (atomstr, parameter.id, parameter.smirks) )
Here we analyze which parameters are applied to a group of molecules.
This example uses get_molecule_parameterIDs, a simple utility function similar to label_molecules, but intended for use on large datasets. get_molecule_parameterIDs processes a list of molecules using a specified ForceField and returns the parameters that would be assigned to each molecule, grouping the results both by molecule and by parameter.
This example may be useful when adding new parameters, to ensure that they are applied to molecules in a given data set.
This example also highlights the difference between "specific" and "generic" parameters. In the SMIRNOFF format, more than one SMIRKS-based parameter may match a motif in a molecule. Therefore, some parameters have precedence over others. Parameters which are "generic" are listed near the top of a parameter section, while more "specific" parameters occur near the end of the section. During parameter assignment, the furthest-down parameter that matches a motif in the molecule is assigned.
During actual parameter assignment, if the toolkit is unable to assign SMIRKS-based parameters to any part of a molecule, an UnassignedValenceParameterException will be raised.
In [2]:
from openforcefield.utils import get_data_file_path
from openforcefield.topology import Molecule
from openforcefield.typing.engines.smirnoff import ForceField
from openforcefield.utils.structure import get_molecule_parameterIDs
# The set of molecules that we want to parametrize in SMILES format.
molecule_smiles = [
'[H][C@@]1([C@](OC1([H])[H])([H])O[H])O[H]',
'[H]C1(C(C(C1([H])O[H])([H])[H])([H])O[H])[H]',
'[H]C1(C(C(O1)([H])[H])([H])O[H])[H]',
'[H][C@@]1(C([C@@](OC1([H])[H])([H])O[H])([H])[H])O[H]',
'[H][C@@]1(C(O[C@](O1)([H])O[H])([H])[H])O[H]',
'[H][C@]1(C([C@](O1)([H])O[H])([H])[H])O[H]',
'[H][C@]1(C(OC(O1)([H])[H])([H])[H])O[H]',
'[H][C@@]1([C@](C(C(C1([H])[H])([H])[H])([H])[H])([H])O[H])O[H]',
'[H][C@]1([C@@](C(OC(C(C1([H])[H])([H])[H])([H])[H])([H])[H])([H])O[H])O[H]',
'[H][C@@]1([C@](OC(C(C1([H])[H])([H])[H])([H])[H])([H])O[H])O[H]'
]
# Create the Molecule objects to parametrize.
mols = [Molecule.from_smiles(smiles) for smiles in molecule_smiles]
# Let's label using the smirnoff99Frosst force field
forcefield = ForceField('openff-1.0.0.offxml')
# This utility function creates dictionaries describing parameter assignment,
# grouped both by molecule and by parameter
parameters_by_molecule, parameters_by_ID = get_molecule_parameterIDs(mols, forcefield)
print("Molecules with parameter IDs:")
# parameters_by_ID is a dictionary where keys are parameter IDs and values
# are molecules in which each parameter ID occurs.
pids_ordered = sorted(list(parameters_by_ID.keys()))
for pid in pids_ordered:
print(pid)
for ids in parameters_by_ID[pid]:
print('\t', ids)