Skip to content

OFF-EP 0010 — Clarify ProperTorsion implementation of idivf and dihedral calculation

Status: Submitted

Authors: Lily Wang

Acceptance criteria: Unanimity (4 approving reviews) or partial support (2 approvals and 2 week period with no reviews requesting changes)[https://openforcefield.atlassian.net/wiki/spaces/MEET/pages/2638774273/09-05-23+SMIRNOFF+Committee+Meeting]

Stakeholders: John Chodera, David Mobley, Jeff Wagner, Lily Wang, Matt Thompson

Created: 2024-02-07

Discussion: Issue #59, Issue #60

Implementation: openforcefield/standards

Abstract

This change clarifies the implementation of SMIRKS matching, dihedral calculation, and idivf="auto" for ProperTorsion parameters.

Motivation and Scope

A ProperTorsion is defined between a connected quartet of atoms i-j-k-l. The dihedral angle is calculated between the two planes defined by i-j-k and j-k-l. The directions of these planes, and the resulting sign of the dihedral, depend on how they are defined; however, it is currently unclear what standard to follow. While both directions yield the same result for symmetric torsions where the phase is 0 or pi, the choice of direction is important for asymmetric torsions.

In addition, a SMIRKS pattern that can match a particular bonded quartet in either i-j-k-l or l-k-j-i order is ambiguous, and the SMIRNOFF specification does not guarantee that the match will be performed in any predetermined or deterministic order. As this may potentially lead to undesired results, this proposal adds a note highlighting this fact.

Finally, the effect of the idivf="auto" parameter on the ProperTorsion potential is outlined in words that can be interpreted ambiguously.

All of these ambiguities may cause confusion to anybody implementing the SMIRNOFF spec.

Usage and Impact

This change defines the idivf parameter in line with the philosophy used in AMBER force fields. The OpenFF implementation has not hitherto used or implemented idivf="auto" in its ProperTorsion parameters. As such there should be no practical impact on existing workflows and simulations. Other force fields and software that have implemented and interpreted the idivf="auto" parameter in ways that do not align with our definition will need to be updated accordingly.

The proposed clarification of input vector order and symmetric SMIRKS matching reflects existing implementations and simulations in OpenFF and OpenMM. Software that converts systems out into other formats may need to adjust the input order to ensure the correct sign of the torsion.

Backward compatibility

This proposal does not change behaviour, but rather explicitly defines what is currently implicit in the specification and the implementation in OpenFF infrastructure. Therefore, there should be no backwards compatibility issues.

Detailed description

This proposal adds the following section (bolded) clarifying idivf to the ProperTorsion spec:

For convenience, an optional attribute specifies a torsion multiplicity by which the barrier height (k#) should be divided (idivf#). The final barrier height is calculated as k#/idivf#. idivf can be assigned an integer value (such as "1"), or "auto". If idivf="auto", the following equation is used to determine the idivf value for a torsion applying to four atoms i-j-k-l, where n_j refers to the degree (i.e. number of bonds) of atom j:

** **idivf = (n_j - 1) * (n_k - 1)** **

The default behavior of this idivf can be controlled by the top-level attribute default_idivf (default: "auto") for <ProperTorsions>.

It also adds a section explaining the computation of theta to the ProperTorsion spec:

In the potential function, the angle theta is calculated using input vectors defined by the four atoms of the torsion i-j-k-l.

Dihedral angle figure

Where the vector r_ij is defined as the vector from atom j to atom i: r_ij = x_i - x_j the angle theta should be calculated using the input vectors r_ij, r_kj, and r_kl. These define the planes u_ijk and u_jkl (see figure below, > section A).

Dihedral angle process

The sign of the angle is determined by comparing the r_ji vector to the u_jkl plane (see figure above, section B). If the r_ji vector has an acute angle to the u_jkl vector, the sign is positive; if the angle is obtuse, the sign is negative (section C in figure above).

Pseudocode of the expected implementation is provided below.

``` u_ijk = r_ji x r_jk u_jkl = r_jk x r_lk angle = acos(u_ijk • u_jkl) # returns in domain [0, pi]

rij_to_ujkl = r_ji • u_jkl if rij_to_ujkl < 0: sign = -1 else: sign = 1 theta = sign * angle ```

Note

Angle values close to 0 and π may be susceptible to precision errors in implementations.

The sign of the theta angle is important in cases where the torsion profile is > asymmetric, i.e. where the phase is neither 0 nor pi, for example in the case below.

Dihedral angle torsion profile

And finally adds a note on how ProperTorsion SMIRKS are applied:

Note

A SMIRKS pattern that can match a particular bonded quartet in either i-j-k-l or l-k-j-i order is ambiguous, and the specification cannot guarantee the match will be performed in any predetermined or deterministic order, potentially leading to undesired and undefined results.

Alternatives

Alternative 1: the idivf parameter could be removed

In this alternative scenario, idivf would be implicitly set to 1, as has been the case in current SMIRNOFF force fields. Torsion parameters would be explicitly enumerated to only apply to one specific multiplicity, and the k barrier would be appropriately fit to the required scale.

The approach of keeping and implementing idivf="auto" was chosen to enable scientific experiments to investigate whether torsion multiplicity could be accounted for using idivf, reducing the number of necessary parameters to fit.

Alternative 2: disallow asymmetric torsions

In this alternative scenario, proper torsion parameters with asymmetric profiles (e.g. with phases outside 0 or pi) would be explicitly disallowed in the SMIRNOFF spec and OpenFF infrastructure. An error would be raised on reading these. This would render concerns about dihedral sign and the impact of atom order superfluous.

The approach of keeping asymmetric torsions was chosen to: * enable reading our existing Sage 2.0 force field, which was published with asymmetric torsions * allow asymmetric torsions for future research

Discussion

This template is based upon the numpy NEP template and the conda-forge CFEP template.

All OFF-EPs are explicitly CC0 1.0 Universal.