smirnoff99Frosst
smirnoff99Frosst copied to clipboard
openforcefield
Differences in how some torsions are modeled in SMIRNOFF vs GAFF for a molecule in FreeSolv
While spot-checking some gas-phase simulations of molecules from FreeSolv parameterized using SMIRNOFF, one molecule that caught our eye in particular was "aldicarb" (ID within FreeSolv: 'mobley_5200358', SMILES: 'CC(C)(/C=N\OC(=O)NC)SC').
Upon further inspection, it appears that there are a few torsions in this molecule that are modeled very differently in GAFF vs SMIRNOFF, but I don't know enough to tell whether the differences are intentional (e.g. because smirnoff fixed issues that were present in gaff), or if any of these require further attention.
Another thing I noticed here was that a lot of the torsion terms in smirnoff99Frosst.offxml have force constants of exactly zero (there are 38 instances of k1="0.000" in the PeriodicTorsionForce definition).
(Code to reproduce this figure: https://gist.github.com/maxentile/6d2610fc72c077d6e605f256f916881f [requires openeye license])
(@davidlmobley suggested to migrate this discussion from Slack : tagging participants from that discussion @jchodera @cbayly13 @mrshirts )
Update: To get a sense of whether this reflects something corner-case-y about aldicarb, I looped this script over the rest of FreeSolv. This reported over 400 molecules where at least one torsion appears to be modeled differently between GAFF and SMIRNOFF.
Here's a random sample of the analogous plots for some of these molecules: torsion-comparison-sample.zip
I was surprised to find so many differences because the SMIRNOFF paper shows that SMIRNOFF reproduces GAFF hydration free energies very closely for virtually all of the FreeSolv set. Maybe these torsion differences shouldn't have been surprising though -- I haven't yet checked how many of these differences are due to human errors in parm@Frosst / GAFF that were corrected in SMIRNOFF, e.g. the torsional energy differences described in section 5.1.1. of the SMIRNOFF paper (which would explain discrepancies for torsions containing atom types H1, H2, or H3).
If this is worth looking into further, I can continue next week.
We expect many differences in torsions; most of the hydration free energies here are not very sensitive to torsions since the molecules are in many cases small and fragment like. If you want to make sure we have access to your scripts it's worth following up on at some point, but not super high priority. Some will likely be cases where GAFF is wrong and some the other way around. And some, both may be wrong. Discrepancies may be compounds to get @chayastern to prioritize for QM torsion scans. Also this may overlap with an energy minimization project we're doing internally (cc @bannanc and @vtlim)
And if i had to guess how many torsions are modeled differently (as in, different barrier heights) in smirnoff99frosst vs GAFF, I would guess that it would be ALL of them except those with barriers of zero in both.
Thanks @davidlmobley for tagging me. I'm slowly catching up on messages from last week. I added some C~S bond parameters from GAFF2, but I just double checked and those were only for double and aromatic bonds (documentation in issue #41). Otherwise I'm fairly sure all of these parameters came from parm99/parm@frosst.
It would be interesting to check what these parameters were in parm99/parm@frosst. That would at least tell us if the difference is "intensional" that is a difference we would expect to see or if there was a typo or issue in porting over the atom types. I don't think we have FreeSolv with parm@frosst atom types though so we would have to find similar functional groups in DrugBank or Zinc in order to do that comparison.
I honestly think the best procedure going forward here would be to have @ChayaSt do some torsion drives on these and head towards fitting them.
I think its worth having chemistry documented that we would like to include in the refitting so this issue is worth including. However, I think as soon as we start fitting parameters that are not just ported from another force field it will be close to time to make a new "SMIRNOFF" that is the point of SMIRNOFF99Frosst is to at least closely replicate parm99/parm@frosst in the SMIRNOFF format.
@yudongqiu have you looked at this in the context of your current refitting?