Torsions with 90/270 degree phases are sometimes mis-applied

[wenyan4work]

Dear OpenFF developers,

When we are using this wonderful package, we found there are three propertorsion parameters in the sage offxml file somewhat mysterious:

<Proper smirks="[#16X2,#16X1-1,#16X3+1:1]-[#6X3:2]-[#6X4:3]-[#7X4,#7X3:4]" periodicity1="4" periodicity2="3" periodicity3="2" periodicity4="2" periodicity5="1" phase1="0.0 * degree" phase2="0.0 * degree" phase3="0.0 * degree" phase4="270.0 * degree" phase5="90.0 * degree" id="t25" k1="0.1024566756215 * mole**-1 * kilocalorie" k2="0.1025622546566 * mole**-1 * kilocalorie" k3="-0.7270841288504 * mole**-1 * kilocalorie" k4="0.2720289395625 * mole**-1 * kilocalorie" k5="0.4515189222906 * mole**-1 * kilocalorie" idivf1="1.0" idivf2="1.0" idivf3="1.0" idivf4="1.0" idivf5="1.0"></Proper>
<Proper smirks="[#16X2,#16X1-1,#16X3+1:1]-[#6X3:2]-[#6X4:3]-[#7X3$(*-[#6X3,#6X2]):4]" periodicity1="4" periodicity2="3" periodicity3="2" periodicity4="2" periodicity5="1" periodicity6="1" phase1="270.0 * degree" phase2="0.0 * degree" phase3="180.0 * degree" phase4="270.0 * degree" phase5="270.0 * degree" phase6="0.0 * degree" id="t26" k1="0.3613352840084 * mole**-1 * kilocalorie" k2="0.7094107693902 * mole**-1 * kilocalorie" k3="0.2966149422494 * mole**-1 * kilocalorie" k4="0.3825271697743 * mole**-1 * kilocalorie" k5="0.01001366302186 * mole**-1 * kilocalorie" k6="-0.05061714024213 * mole**-1 * kilocalorie" idivf1="1.0" idivf2="1.0" idivf3="1.0" idivf4="1.0" idivf5="1.0" idivf6="1.0"></Proper>
<Proper smirks="[#6X4:1]-[#16X4,#16X3+0:2]-[#7X3:3]-[#6X3:4]" periodicity1="3" periodicity2="2" phase1="90.0 * degree" phase2="0.0 * degree" id="t151" k1="-0.6373653602282 * mole**-1 * kilocalorie" k2="0.2973949413807 * mole**-1 * kilocalorie" idivf1="1.0" idivf2="1.0"></Proper>

More specifically, they are:

1. id="t25" phase4="270.0 * degree" phase5="90.0 * degree"
2. id="t26" phase1="270.0 * degree" phase4="270.0 * degree" phase5="270.0 * degree"
3. id="t151" phase1="90.0 * degree"

These phase values transform the cos series from an even function to an odd function. Therefore, when a dihedral angle is calculated from the atom indices (ordered) as 1-2-3-4 or 4-3-2-1, we get different results.

The label_molecule function in openff toolkit seems to always return atom indices according to their indices in the Molecule data, therefore, the order of the returned 4 atom indices does not always match the order in the defined smirks pattern.
For example, for t25, the smirks pattern defined the atom sequence as S-C-C-N, but in a molecule the N atom may appear before the S atom, so in the returned (sorted) atom indices, this group of torsion may be listed as N-C-C-S, and cause some discrepancies for the torsion energy calculations.

Did I misunderstand the 'label_molecule' function or the role of the 'phase' parameter?

Thank you,

[wenyan4work]

I found these phase angles come from this file: http://www.ccl.net/cca/data/parm_at_Frosst/parm_Frosst.frcmod and has been used in all versions of openff offxml files

Also, in 0.10.6 (not in 0.11.0), the two permutations of the propertorsion atom indices look weird. Is this a bug or a feature?
https://github.com/openforcefield/openff-toolkit/blob/37d6806bbb9d0d8aacb7221891b5d86dcdee014b/openff/toolkit/topology/topology.py#L146

Should this be 3-2-1-0, instead of 3-1-2-0?

[lilyminium]

Should this be 3-2-1-0, instead of 3-1-2-0?

Just to comment on this -- that's right, it should be 3, 2, 1, 0 @wenyan4work. It was a minor bug, fixed in openforcefield/openff-toolkit#1185. However, I'm not sure ValenceDict.index_of is/was actually used anywhere in the code.

[wenyan4work]

"[#16X2,#16X1-1,#16X3+1:1]-[#6X3:2]-[#6X4:3]-[#7X4,#7X3:4]"

Thank you, @lilyminium !

Could you also comment on the issue of the order of atoms in the result returned by 'label_molecule'? For example in the parameter t25, the order defined by the smirks pattern is S-C-C-N, but in our molecule since the N has an atom index smaller than S, the parsed result becomes N-C-C-S.

[j-wags]

This is a really great question @wenyan4work.

I could be wrong, but I don't really see a physical justification why general dihedral parameters should have a "handedness" - That is, I'd think that all of our dihedral parameters should be "even". Chemically, if we intended parameters to have a handedness, I'd also expect to see stereochemistry in the SMIRKS (or for the SMIRKS to be constructed such that they guarantee some symmetry, so that the term is applied many times across the central bond and the handedness cancels out)

In addition to phase, each ProperTorsion term also has a periodicity (basically a frequency). I was thinking that maybe these were constructed with a periodicity that interacts with the phase to remove the handedness - For the torsions you've identified they are:

• id="t25"
  • periodicity4="2" phase4="270.0 * degree"
  • periodicity5="1" phase5="90.0 * degree"
• id="t26"
  • periodicity1="4" phase1="270.0 * degree"
  • periodicity4="2" phase4="270.0 * degree"
  • periodicity5="1" phase5="270.0 * degree"
• id="t151"
  • periodicity1="3" phase1="90.0 * degree"

Plotting some of these as cos((periodicity*x) - phase), some of them look even, but others are indeed odd.

The other weird thing about these parameters is that they have multiple terms with the same periodicities - t25 and 26 have periodicities of 4,3,2,2,1 and 4,3,2,2,1,1 respectively. I don't understand the purpose of the repeated periodicities.

I'm looking through other information sources and will compile information here - This isn't the first time this has come up but I'm not sure that the question was resolved before. Most conversations see that the 90 and 270 values were in parm@Frosst, and everyone assumes they must have a good reason, and the discussion ends there. But maybe it was a mistake all along.

[j-wags]

Slack DM thread from March 2021:

Christopher Bayly
13:39
We are puzzled by torsions having more than one parameter of the same periodicity, e.g. t25 in 1.3.0 which has two periodicity 2 torsions (they have different phase angles):

It was unexpected by us so we are fixing our related bug but it seems really really odd to me (physics-wise). Without using too much of your time can you tell me when we started doing this?

Jeffrey Wagner
14:20
Looking quickly, it looks like the doubly-occupied periodicity of t25 started in S99F-1.0.7 [1]. But even before then (in version 1.0.6 and earlier), t26 and t27 has the same thing going on
[1] https://github.com/openforcefield/smirnoff99Frosst/blob/master/smirnoff99frosst/offxml/smirnoff99Frosst-1.0.7.offxml#L167
[2] https://github.com/openforcefield/smirnoff99Frosst/blob/master/smirnoff99frosst/offxml/smirnoff99Frosst-1.0.6.offxml#L169-L171

Christopher Bayly
14:22
Huh... did I do that, way back in the day? It doesn't look like the kind of thing I would do...
14:22
Thanks for looking into it.

[j-wags]

cc martimunicoy/peleffy#54