...
While OpenFF has yet to move to a full neural network force field in the framework of Espaloma, it may be useful to use Espaloma as a reference, and we may be able to use Espaloma to determine areas where OpenFF parameters need improvement. For example, there may be cases where OpenFF uses one parameter to encode a particular chemistry, that Espaloma splits into many different values. Here, Trevor Gokey’s work on automated parameter generation could come in handy for partitioning espaloma data. If assigned parameter values are significantly different between Espaloma and OpenFF, that would also be worth exploring.
| Table of Contents |
|---|
Data
This PDF contains histograms comparing all of the espaloma-assigned values for a given Sage parameter to the value of the Sage force constant.
| View file | ||
|---|---|---|
|
However, a much more useful way of looking at this data can be found in this repository. After installing the corresponding conda environment, running make will build the necessary data files, and then running either python board.py or python twod.py will allow you to to view the data interactively in a browser. Below is an example of the interface taken from the README.
...
Espaloma Benchmark
The figure below shows a benchmark comparing Sage 2.1.0, the Sage 2.1.0 force field with many force constant values taken from the average Espaloma values, and Espaloma itself. Espaloma performs very well.
...
Experiments
Big Torsion Deviations
As a first attempt, I labeled a data set with both Espaloma and Sage 2.1.0 and compared the values they assigned. Two of the torsions, shown below, had deviations between the Sage force constant and the average Espaloma value of more than 10 kcal/mol. These correspond to torsion IDs t129 and t140, respectively.
...
For these torsions, I replaced the Sage value with the average value from Esplaoma in the force field, and ran benchmarks on the OpenFF Industry Benchmark Season 1 v1.0 data set, yielding the plots below. I didn't expect to see much difference from such a small change of only two parameters, but it's encouraging that it didn't ruin anything, at least. The eps-tors-10 results might even be very slightly better, as desired.
...
All Parameters
With these results in hand, I next repeated the process but replacing every Sage parameter with the corresponding average parameter from Espaloma. This probably isn’t the best approach because many of the distributions look like the one shown below: there are multiple clusters of Espaloma-assigned values, and the Sage value is out in the middle. These may be good candidates for parameters that need to be split in Sage.
...
As shown below, the results are more different from the esp-tors-10 results, as expected. And positively, esp-full appears to perform a bit better by all three metrics. This is without any re-fitting, so Espaloma’s average parameters for our SMIRKS patterns perform slightly better than our re-fit Sage 2.1.0 values.
...
Possible Splits
b84
...
The figure above shows the distribution of Espaloma parameters for b84. Whereas the Espaloma average is at 741.5, the Sage value is dragged down to 719.6 by a small cluster of pattern matches near 700. Zooming in on this cluster shows that nearly all of the carbons are bound to nitrogen, with 3/104 bound to oxygen instead. This suggests that this parameter could potentially be refined by adding a more specific parameter like [#7X3]-[#6X4:1]-[#1:2] or [#7X3,#8X2]-[#6X4:1]-[#1:2]. Examples of these molecules are shown in the figures below.
...
b5
...
For b5, the long tail of values below the Sage value is due to C-C bonds between two fused aromatic rings, as shown in the figure below. It would probably be difficult to split b5 enough times to replicate all of the different Espaloma values along this tail, but it might be reasonable to split it at least once to separate the main body of Espaloma values centered at 843 from the tail.
...
b1
...
Initially, I thought these were mostly CX4-CX4 bonds in rings. Many of them are actually in three- or four-member rings, but there are some more typical C-C bonds too. There are actually 396 of these, which was a bit too much to crawl through by hand. A random selection of 25 are shown below. As a result of the variety, I’m not sure this one is a great candidate for splitting. Maybe we’re looking for distributions with two (or more) distinct peaks rather than ones with long tails like b1 and b5.
...
b2
...
As shown in the figure below, all of these outliers are due to C-C bonds where one or both C is attached to an O, usually a carbonyl or carboxyl group.
...
b7
...
b16
...
b4
...
b21
...
b6
...
b10
...
b9
...
b13
...
b59
...
b3
...
b19
...
b11
...
b52
...
b88
...
b17
...
b14
...
b20
...
b36
...
The results above were actually still relying on the code that identified the large deviations from Sage parameters, so not all parameters were replaced. Additionally, that code only looked at force constants, and only the first force constant for torsions. In contrast, the graph below shows the esp-full-full (really full) results, where all of the values have been replaced for all force constants and equilibrium bonds and angles. Clearly these look much worse for the Espaloma values within the SMIRNOFF format. This makes me a bit suspicious of the torsions especially because those values look the most different from Sage at first glance, at least compared to the bonds and angles, which are more recognizable.