Is your feature request related to a problem? Please describe.
@ChayaSt's work showed that we may be able to produce higher-quality proper torsion parameters by taking into account the "partial bond order" of the torsion's central bond. We now have the machinery to compute AM1-Wiberg partial bond orders for entire molecules using either OpenEyeToolkitWrapper's or AmberToolsToolkitWrapper's assign_wiberg_bond_orders methods. The thought is that, if some simple electron population analysis shows that a bond order is 1.53, maybe torsions about that bond can be described well by interpolating to 53% of the way between the single and double bond k values.

A similar area of the spec has been proposed for interpolation of parameters in the Bonds section. These keywords are currently unsupported, but provide a template we can follow for torsion interpolation. If needed, we could implement Bond interpolation as a technical exercise, but that functionality is far lower-value than proper torsion interpolation.

Describe the solution you'd like

The extension of the ProperTorsions tag in the SMIRNOFF XML spec, adding two new supported ProperTorsions header attributes:

• fractional_bondorder_method, with default value AM1-Wiberg
• fractional_bondorder_interpolation, with default value linear

The extension of the ProperTorsion tag in the SMIRNOFF XML spec to support EITHER

• (currently supported) a set of non-interpolation-based torsion parameters, with attributes (for nonzero, consecutive integer values of X):
  • kX (required)
  • phaseX (required)
  • periodicityX (required)
  • idivfX (optional)
• (to be added) a set of a set of interpolation-based torsion parameters, with attributes (for nonzero, consecutive integer values of X and nonzero, consecutive integer bond orders B):
  • kX_bondorderB (required)
  • phaseX (required)
  • periodicityX (required)
  • idivfX (optional)

Note: I have no strong leaning as to whether, for a single parameter, kX_bondorderB should support a different range of B values for different Xs. We can let ease of implementation guide this.

Note: I have no strong leaning as to whether B values should be required to be integer, however it's quite clear that this will make implementation much easier, so I think it's the way to go for our first attempt.

Note: We may consider what to do if a bond order falls outside the range of B values. We could solve this by adding an extrapolate attribute to the ProperTorsions tag, which could default to False (in which case, a bond order of 0.91 would get the k values from kX_bondorder1), but the FF could set it to True (in which case, a bond order of 0.91 would be linearly extrapolated from the kX_bondorder1 and kX_bondorder2 values)

Implementation details

We will want to add a keyword argument to ForceField.create_openmm_system called partial_bond_orders_from_molecules=None, but which could be set to [iter of OFFMol] (similar in implementation to the existing charge_from_molecules keyword argument), which lets users override our automatically-generated partial bond orders if they wish.

In ProperTorsionHandler.create_force, adding logic to run assign_fractional_bond_orders only if necessary, since a molecule might not be assigned any fractional bond order-dependent parameters, or the user may have specified partial_bond_orders_from_molecules, and we wouldn't want to waste time computing them in that case.

In the class definition for ProperTorsionHandler.ProperTorsionType, we'll need to either:

• define a new subclass of ParameterAttribute, like DoubleIndexedParameterAttribute.
• somehow change or subclass IndexedParameterAttribute to be able to "wrap" and add a SECOND indexing system to an existing IndexedParameterAttribute.

I'm happy to discuss further details or consider changes to this specification. Since this will be our first version of support for torsion interpolation, we don't have any existing behavior to maintain!