Improving Small Molecule Force Fields by Identifying and Characterizing Small Molecules with Inconsistent Parameters

2.1. Molecule pairs were flagged as similar or different based on TFD and TanimotoCombo

We used TanimotoCombo and TFD to identify molecules with dissimilar geometries to seek molecules with parameter inconsistencies. We assign a “difference flag” to a molecule pair (in a “molecule pair”, the comparison is made across force fields) when it yields a TFD value over 0.20 and a TanimotoCombo value over 0.50. These pairs visually exhibit different minimized geometries that may be indicative of parameterization differences. Out of 26,984,560 possible molecule pairs involving any pair of force fields, the combination of the SMIRNOFF99Frosst and GAFF2 force fields yielded the largest number of difference flags (305,582, Table 1). This indicates that these force fields are quite different. In contrast, the combination of MMFF94 and MMFF94S yielded the smallest number of difference flags at 10,048 difference flags, indicating that these two force fields are the most similar among those being compared. These numbers are sensible given the history of these force fields – GAFF2 has undergone considerable recent reparameterization [39], and SMIRNOFF99Frosst inherits parameters from parm@Frosst [1], a sibling force field of GAFF, while reducing the number of parameters with an entirely different form of chemical perception [26, 27]. In contrast, MMFF94 and MMFF94S are identical aside from their treatment of some nitrogen atoms [15]. Consequently their optimized conformers should be rather similar, as reflected in our scores. Thus, these results match what would be expected from the parameterization history of these force fields.

We also label molecule pairs with highly similar geometries. To do this, we assign “similarity fags” to molecule pairs that yielded TFD values under 0.18, indicative of similar geometries (Table 2). In order to visualize the number of molecule pairs with each fag, we plotTFD versus TanimotoCombo for all molecule pairs in Figure 1. We highlight regions fagged as similar and different along with regions outside the interest of this analysis. Figure 1 likewise shows that the vast majority of molecule pairs were rated similar by both TFD and TanimotoCombo.

2.2. Sets of molecules were created based on their similarity and difference flags

We then sort the molecules into sets of interest by their patterns of difference and similarity fags. As molecule pairs were formed from a set of five conformers, each resulting from optimization with a different force field, each molecule results in ten different molecule pairs which can be assigned either a difference or similarity fag. All molecules that yielded five or more difference fags out often were added to the set named “FivePlus.” We also categorized molecules of particular interest for each force field. For each force field, we identified molecules in which two conditions held: (1) all molecule pairs involving that force field were fagged as different, and (2) the molecule pairs not including that force field were fagged as similar. Accordingly, molecules in these sets must result in four difference fags and six similarity fags; molecules in these sets can not also be in the FivePlus set. This allows us to highlight molecules which were treated differently by only one force field, potentially indicating problems in the force field’s parameters for the represented chemistries of the molecule. We called this set the “Individually Different” set for that force field. For example, the molecules identified in this scheme for SMIRNOFF99Frosst were added to the “Individually Different SMIRNOFF” (ID_SMIRNOFF) set. This latter analysis is probably most relevant to the SMIRNOFF force field, as GAFF/GAFF2 and MMFF94/MMFF94S come in families which would reduce the number of cases meeting these criteria if intra-family similarity is high – specifically, if both family members treat a molecule consistently, it will not be fagged as “individually different” for that force field.

Our results after categorizing put 111,162 molecules into the FivePlus and 93,859 molecules in the ID_SMIRNOFF set out of a total of 2,698,456 molecules. The ID_SMIRNOFF set was the largest of the individually different force field sets, as is displayed in Table 3. As noted, we had some expectation SMIRNOFF might be relatively distinct from the other force fields considered.

Here, we focused on identifying molecules with significant geometric differences between force fields, and our sets were constructed to help identify these molecules, but other factors might also be important to examine in future work. For example, if different force fields lead to similar optimized geometries, that does not necessarily mean those force fields are similar. To examine whether energetics of the different force fields are similar, we would need to study the relative energetics of conformers of different molecules in different force fields, which is not something within the scope of this work as it would require multiple conformers per molecule. However, relative energetics have been examined in a separate study [24]. Here, then, we focus on identifying geometric differences which likely imply force field differences, though geometric similarities do not necessarily imply force field similarities.

2.3. Certain functional groups are more likely to appear in molecules with geometric differences

2.4. This work has been used to improve training datasets for the OpenFF Parsley series

In the present work, discrepancies between optimized geometries from different force fields highlight potential issues, but we have no ground truth or point of reference for sorting out which geometries are correct and which are not. This data simply helps us select molecules/chemistries which may be informative, and priloritize them for further study. Particularly, one might generate optimized geometries for these same molecules with QM calculations and then use these to help assess which force fields produce the best results, or use these in force field training sets to improve force field quality.

Indeed, informative molecules from the present study are being used for precisely that purpose. Particularly, a subset of the FivePlus set was used as the basis for the “coverage” set used for the first OpenFF Parsley release, OpenFF 1.0 [29]. A larger portion was used in benchmarking OpenFF 1.0. Then, for OpenFF 1.2, training data was completely redesigned, in part drawing from what was called the “eMolecules Discrepancies Set” [22, 25], corresponding to the first portion of the FivePlus set generated here. This training data redesign resulted in improved performance on a variety of benchmarks [21, 24]. The relevant optimized geometries are freely available in QCArchive [34] as part of the OpenFF 1.2 training and benchmarking datasets.

While subsequent OpenFF work building on the data generated here is not formally part of this study, it does appear that molecules identified as potentially informative by this approach do serve well as input for QM calculations and force field training, at least when coupled with additional data selection and curation steps.

Multiple force fields were used to minimize conformers

We created input files for multiple force fields from a filtered eMolecules set (filtering described in Section 3.3). We generated molecules from the SMILES strings as in eMolecules, adding explicit hydrogens and assigning default protonation states using the OpenEye toolkits. We did not enumerate protonation states or tautomers, and no significant effort was invested in selecting protonation states; we simply took the default states provided by the toolkit. We do not see this as a major limitation in a force field comparison since the resulting approach tests the force fields thoroughly on the molecules and protonation states used, even if that protonation state or tautomer will not be the most populated at neutral pH in solution.

Following construction of initial molecules, initial conformers were generated with OpenEye’s Omega, then partial charges were assigned to molecules before minimization using the OpenEye implementation of AM1-BCC [19, 20]. The input generation process yields one Tripos MOL2 file to be minimized directly with SMIRNOFF99frosst, MMFF94, and MMFF94S, as well as individual input coordinate and parameter topography files for use by GAFF (1.8) and GAFF2 (2.1). These force fields were chosen because they are widely used, easily available, and compatible with our workflow. Other force fields were either incompatible with our toolchain without substantial additional work, or were commercial and proprietary. For example, comparisons with CGenFF [36, 37], OPLS-AA [23], or the Schrödinger OPLS series [17, 30] would be of considerable interest, but these require substantially different toolchains, and the most recent Schrödinger force fields are also proprietary and require paying for a license.

We minimized each molecule using the parameters from each of the five aforementioned force fields, making sure to start all five minimizations from the same conformer. Minimizations with force fields other than MMFF were performed with OpenMM 7.0.1 [4] using the L-BFGS algorithm [28] with an energy tolerance of 5.0e-9 kJ/mol and a maximum of 1500 iterations. MMFF minimizations were performed with OpenEye’s Szybki Toolkit [35, 42]. Sample run files can be found in the Supporting Information. Molecules that did not successfully result in five minimized structures (one from each force field), were removed from analysis. For each molecule with five minimized structures, pairwise comparisons yielded a total of ten molecule pairs for geometric evaluation. We call these pairs of minimized conformers generated by different force fields “molecule pairs.”

Molecule pairs were assessed using Torsion Fingerprint Deviation and TanimotoCombo

We then assessed each molecule pair for geometric differences. Molecule pairs were evaluated using two distinct measurements: Torsion Fingerprint Deviation (TFD) and TanimotoCombo.

TFD is a method of measuring geometric differences between two conformers of the same molecule based on torsion angles. The TFD score between two structures represents a weighted sum of torsional differences as defined by Schulz-Gasch et al. [31]. Torsions central to the molecule are given more weight than torsions on the periphery of the molecule. Similarly with RMSD, geometric similarity is inversely correlated with TFD score. TFD scores range from 0 to 1, with 0 being most similar and 1 being most different. The authors of TFD consider scores over 0.2 to represent significantly different geometries. In contrast to RMSD, TFD is bounded and less sensitive to molecular size, making it particularly helpful here.

TanimotoCombo, from OpenEye Scientific, is a normalized method of measuring geometric similarity between molecules. It is the sum of ShapeTanimoto, a measure of overall spatial overlap between two molecules, and ColorTanimoto, a measure of spatial overlap of specific functional groups between two molecules, both of which are also metrics from OpenEye. TanimotoCombo values between two conformers range between 0 and 2 (it is the sum of two values each running from 0 to 1), with 2 being the most similar and 0 being the most different.

By visual inspection, we determined that TanimotoCombo is useful for recognizing cases where geometric differences are caused by particularly flexible moieties, such as single bond rotations in an alkyl chain. These differences can often be attributed to minor differences between force fields leading to flexible bond rotations, not to larger differences in force fields that result in more substantial geometric differences. Thus, here, we find that TanimotoCombo alone does not serve to help us isolate geometry differences that are likely due to substantial force field differences; instead, low TanimotoCombo values can result from simple bond rotations that result from molecules energy minimizing to different local minima that we do not consider particularly interesting by visual inspection. However, TanimotoCombo in conjunction with TFD can be used to identify geometric differences that suggest underlying inconsistencies in parameterization.

Molecule pairs were flagged as similar or different based on TFD and TanimotoCombo

We identified molecule pairs displaying parameterization differences which led to different geometries using TFD and TanimotoCombo. TFD is sensitive to ring deformations, torsional differences, and atom planarity changes, which makes it useful for recognizing differences in parameterization. TanimotoCombo, with greater sensitivity to coordinate differences caused by conformational flexibility in a molecule, is more useful for removing cases that are less likely to be caused by parameterization differences, such as different rotameric states.

We chose cutoffs to identify molecule pairs displaying parameterization differences (flagged “different”) and pairs displaying no parameterization differences (flagged “similar”). TFD values below 0.20 are believed to be pharmacologically similar [31], so we chose a TFD value greater than 0.20 to label molecule pairs as different. After visual inspection of a variety of molecules, we observed that molecule pairs with a TanimotoCombo under 0.5 typically had changes due to single bond rotations. Because such bond rotations can arise from a variety of reasons aside from substantial differences in parameterization, we did not wish to focus on such cases. Thus, molecule pairs with a TFD value greater than 0.20 as well as a TanimotoCombo value greater than 0.50 were flagged as different – allowing us to focus on cases with substantial torsional differences which were not simply due to rotations around highly flexible bonds. We used a substantial amount of manual inspection of these thresholds to help us make these choices. As a result of these choices, any pair of molecules with a TFD value of 0.18 or less was assigned a similarity flag, as it will display geometrically similar structures. We left a small buffer region between 0.18 and 0.2 when defining similarity flags in order to avoid an extreme sensitivity to small changes around the 0.20 cutoff.

Molecule pairs that yielded very high TFD or very low TanimotoCombo values were also determined to often be uninformative. Tagging these molecule pairs as “different” would be unhelpful because the differences are not due to substantial changes in force field parameters. Most molecule pairs in this category displayed cases of what might be called “conformer chirality” – where an achiral molecule was minimized to two similar but non-superimposable structures – essentially, collapsing down to two minima which are equivalent to the force field but not geometrically identical. A number of molecule pairs, most with small flexible ring systems, yielded TFD values greater than 1 (which should not be possible). While this behavior was unexpected, we continued to use RDKit’s implementation without modifications. Considering these results, we removed molecule pairs with a TFD greater than 0.60 or a TanimotoCombo less than 0.25.

3.1. We created and characterized sets of interest

Molecules can be sorted into sets of interest by considering the combinations of their difference and similarity flags. A single molecule in this pipeline is associated with five minimized structures. Pairwise combinations of these structures will yield ten molecule pairs and thus up to ten flags. Molecules that yielded a large number of difference flags, regardless of the force fields of origin, are of particular interest for force field parameterization. Specifically, we set aside molecules with five or more difference flags for further analysis, we call this our FivePlus set.

The other sets of interest are based on the origin of the difference flags with the goal of identifying molecules which behave differently with one force field than all the others. For a molecule to be considered different with that one force field, all four molecule pairs involving that force field should be flagged as different, and all other molecule pairs need to be flagged as similar. We call these Individually Different Sets for each force field, i.e. for SMIRNOFF we create the SMIRNOFF Individually Different set labeled by ID_SMIRNOFF. A molecule in the ID_SMIRNOFF set would have 4 difference flags, one for each pair involving SMIRNOFF, and six similarity flags for all other force field combinations.

3.2. Sets of interest were analyzed by the frequencies of the functional groups

Identifying functional groups which are more prevalent in our sets of interest could be informative for future force field parameterization. To this end, we used Checkmol [9] to describe the combination of functional groups in each molecule. When given a molecule, Checkmol provides a list of descriptors for the functional groups it contains. For each descriptor, we count the number of affiliated molecules in each set of interest as well as in the entire molecule set. From there, we can determine the most over-represented descriptors in each set of interest. We only considered descriptors and descriptor pairs that appeared at least twice in our full molecule set.

We compute the over-representation factor describing how over-represented a particular descriptor is in a given set by dividing the frequency of the descriptor in the set by the frequency of the descriptor in the full molecule set. Mathematically, we can write

where N_A,set is the number of molecules containing descriptor A in a particular set, N_mols,set is the number of molecules in that particular set, N_A,total is the number of molecules in total with descriptor A, and N_mols,total is the number of molecules in total.

Force field behavior could change with combinations of functional groups, and thus we repeated this calculation with pairs of Checkmol descriptors. We can apply Equation 1 to analyze pairs of descriptors by replacing A with A + B to represent molecules containing both descriptors. However, with pairs of descriptors, we are more interested in whether the combination of the descriptors is important. For example, if both descriptors A and B are highly probably in a set of molecules, then finding the combination in that set at a higher frequency is not particularly interesting. Thus, we try to determine if the descriptor pair is more likely to show up in a set of interest than the individual descriptors separately. To that end we calculate an enrichment factor given by

where p_{A + B} denotes the observed frequency (probability) of a molecule with the combined A and B descriptors being found in the set of interest, and p_A and p_B denote the individual frequencies for descriptors A and B in the same set of interest. For example, p_A is given by

A larger enrichment factor indicates that the combination of descriptors A and B is more likely to occur in a set of interest than those descriptors individually. Descriptor pairs with a larger enrichment factor should be considered as important for future parameterization because the combination of functional groups changes a force field’s behavior.

3.3 Molecules were sourced from the eMolecules online database

Approximately 8.1 million molecules were initially sourced from the eMolecules database as SDF files (version obtained in September 2016) [5]. Molecules from this set were then filtered by several criteria. We removed all molecules that contained any metal or metalloid atoms, were over 200 heavy atoms, or had a nonphysical valence (such as a pentavalent carbon atom). Molecules which failed at any step of the process were also removed, i.e. could not be parameterized by one of the force fields. While we minimized all these molecules with each force field, very large molecules are impractical for visual inspection or future QM calculations. Thus, we filtered the molecules for analysis here to remove molecules with more than 25 heavy atoms.