New out-of-plane functional forms for fitting experiments

[EvienCheng]

Description

Adds three new out-of-plane potential forms as SMIRNOFF plugins (HarmonicHeight, HarmonicAngle, LeeKrimm), alongside the existing SMIRNOFF periodic-torsion treatment, giving a choice of functional forms for out-of-plane deformation at trivalent centers.

Co-authored-by: @pavankum , and a final code review and the messy git history sorted by Claude.

Background: out-of-plane functional forms

SMIRNOFF (periodic torsion). Out-of-plane deformation is treated as a periodic torsion over the dihedral $\theta$ between the $(i,j,k)$ and $(j,k,l)$ planes:

Depiction of SMIRNOFF's three improper torsions with the same handedness, out of six possible permutations with the central atom (second atom) fixed. Improper torsions enforce planarity, so the ideal out-of-plane angle is 0°, and the potential must penalize deviations symmetrically. A periodicity of 2 and phase of 180° places the energy minimum at $\theta$ = 0° and ensures symmetric penalties for bending above or below the plane.

$$E_{OOP} = k (1+\cos(\text{periodicity}\cdot\theta - \text{phase}))$$

Alternate functional forms being added here

Harmonic, Wilson angle (HarmonicAngle). Uses the classic Wilson–Decius–Cross out-of-plane angle $\theta$, symmetrized over the three angles at a center:

$$E_{OOP} = \tfrac{1}{2} k (\theta-\theta_0)^2$$

Here $\theta$ is the Wilson out-of-plane angle, defined as the angle between a bond $\mathbf{x}_{41}$ and the plane containing the other two bonds $\mathbf{x}_{42}$ and $\mathbf{x}_{43}$:

$$\sin\omega_1 = \mathbf{e}_{41}\cdot\left(\frac{\mathbf{e}_{42}\times\mathbf{e}_{43}}{\sin\phi_1}\right)$$

where $\phi_1$ is the valence angle between $\mathbf{x}_{42}$ and $\mathbf{x}_{43}$. In practice, the symmetrized average of the three Wilson angles at a given center is used. This is the classic definition from Wilson, Decius & Cross and remains widely used in molecular mechanics force fields such as MMFF94 and MMFF94s, which adopt the form:

$$E_{OOP} = 0.043844,\frac{k^{oop}_{IJK:L}}{2},\chi^2_{ijk;l}$$

where $\chi_{ijk;l}$ is the Wilson angle in degrees. The harmonic approximation around $\theta_0 = 0^\circ$ is appropriate for most sp$^2$ centers near planarity.

Harmonic, pyramid height (HarmonicHeight). Uses the dimensionless pyramid height $h$ (Lee, Palmo & Krimm), the normal distance from the apex atom to the plane of its three neighbors (figure below from the paper above):

$$E_{OOP} = \tfrac{1}{2} k (h - h_0)^2$$

The coordinate $h$ is defined as:

$$h = \hat{v} \cdot \mathbf{e}_{41} = \frac{1}{v},\mathbf{e}_{41}\cdot(\mathbf{e}_{42}\times\mathbf{e}_{43})$$

where $\mathbf{e}_{41}$, $\mathbf{e}_{42}$, $\mathbf{e}_{43}$ are unit bond vectors and $\hat{v}$ is the unit normal to the plane spanned by the three end atoms. The coordinate $h \in [-1, 1]$, with $h=0$ corresponding to a planar structure. Unlike the Wilson angle, well-defined for all admissible geometries with no extra symmetrization needed. $h_0=0$ for planar ($sp^2$) centers; $k$ fit to local curvature.

Lee–Krimm (LeeKrimm). Built on pyramid height $h$, designed to fix both of the above:

$$E_{OOP} = V_2\left(\frac{|h|^t}{1-|h|^s}\right)^{2} + V_4\left(\frac{|h|^t}{1-|h|^s}\right)^{4}$$

The parameters $V_2$ and $V_4$ are the quadratic and quartic spring constants for small and large puckers, respectively, while $t$ controls the curvature of the potential near $h=0$ and $s$ (a positive integer) governs how steeply the potential rises as $h\to\pm 1$; larger values of $s$ make the function behave like a polynomial in $h^2$ near $h=0$.

Test plan

• Unit tests for UreyBradley, HarmonicHeight, LeeKrimm, HarmonicAngle collections/handlers pass
• Verified against methane test geometries

Status

• Ready to go

[pavankum]

(Not urgent!)

Hi @fjclark & @mattwthompson , here is some work from Evien, an undergrad in our lab, whenever you have time can you please review this, thank you!

[mattwthompson]

Thanks @EvienCheng - I will have a look through the code next week

[codecov-commenter]

Codecov Report

❌ Patch coverage is 89.47368% with 16 lines in your changes missing coverage. Please review.
✅ Project coverage is 87.05%. Comparing base (8988f83) to head (75af216).

Files with missing lines	Patch %	Lines
smirnoff_plugins/collections/valence.py	86.55%	16 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #110      +/-   ##
==========================================
+ Coverage   86.59%   87.05%   +0.45%     
==========================================
  Files           9        9              
  Lines         694      842     +148     
==========================================
+ Hits          601      733     +132     
- Misses         93      109      +16     

Flag	Coverage Δ
unittests	87.05% <89.47%> (+0.45%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.