Sigma parameters in 1-4 scaling factors are always computing using Lorentz-Berthlot mixing rules

[mattwthompson]

Bug summary

This one goes a little deep, so buckle up. Since Foyer's force field class does not explicitly encode the mixing rule/combination rule, it's somewhat ambiguous how it is ultimately applied. If the force field is meant to encode Lorentz-Berthelot mixing rule, things should work well because that is ParmEd's default and OpenMM assumes that throughout using openmm.NonbondedForce. If using geometric mixing rules, I believe the normal flow is to call Forcefield.apply() on something and then use the setter of the resulting ParmEd Structure, i.e. my_struct.combining_rule == "geometric". Then, ParmEd's writers will generally handle what needs to be handled. For example, writing out to a GROMACS topology file will set comb-rule appropriately (3 is geometric):

;   File foyer.top  was generated
;   By user: mwt (501)
;   On host: Clausius.local
;   At date: Fri. June  5 14:57:26 2021
;
;   This is a standalone topology file
;
;   Created by:
;   ParmEd:       test_foyer.py, VERSION 3.4.1
;   Executable:   test_foyer.py
;   Library dir:  /Users/mwt/miniconda3/envs/openff-system/share/gromacs/top
;   Command line:
;     test_foyer.py
;

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               3               no              1            0.5

This tells GROMACS to use geometric mixing rules when computing the cross-interactions for each atom type. A deep subtlety here is that the cross-interactions of 1-4 pairs have already been computed by ParmEd based on the values given to it from the openmm.System in the pmd.openmm.load_topology() call here. To repeat, openmm.NonbondedForce constructs everything assuming Lorentz-Berthelot mixing rules and encoding geometric mixing rules takes a little more work (https://github.com/openmm/openmm/issues/1859, or @ahy3nz has a deep dive that covers this point: https://ahy3nz.github.io/posts/2019/30/openmm2/). Included in an openmm.System is the exceptions (OpenMM language), also knows as adjusts (ParmEd/Amber? language) or pairs (GROMACS language). So the sigma values in the ParmEd Structure.adjusts list have been pre-computed incorrectly.

I think I've reproduced this below, but the resulting errors are small (0.02 kj/mol/molecule) that I could easily be making a mistake that invalidates my conclusions.

I don't think fixing this elegantly would be possible; the easiest way to stop the bleeding would probably involve adding behavior around here that, if using a geometric mixing rule, goes through each adjust.type.sigma in structure.adjusts and overwrites the existing value with a freshly computed geometric mean. A long-term fix involves a lot of OpenMM hacking that I don't think ParmEd is designed to adjust; I can expand on that if desired but it's no small task. Explicitly storing the combing rule in Foyer's force field object would be necessary for any sustainable solution IMHO.

Code to reproduce the behavior

In:

from typing import TYPE_CHECKING

import foyer  # type: ignore
import mbuild as mb  # type: ignore

if TYPE_CHECKING:
    from parmed import Structure

def _patch_14_with_geometric(struct: "Structure") -> "Structure":
    """Modify sigma values in Structure.adjusts to follow geometric mixing rules"""
    for adj in struct.adjusts:
        sig1 = adj.atom1.sigma
        sig2 = adj.atom2.sigma
        sig_geo = (sig1 * sig2) ** 0.5
        adj.type.sigma = sig_geo

    return struct

oplsaa = foyer.forcefields.load_OPLSAA()
benzene = mb.load("c1ccccc1", smiles=True)

out = oplsaa.apply(benzene)

assert out.combining_rule == "geometric"

# Extract a C-H 1-4 pair from the Structure.adjusts list
for adj in out.adjusts:
    if adj.atom1.name == "C" and adj.atom2.name == "H":
        break

found_14_sigma = adj.type.sigma
sigma_if_lorentz = (adj.atom1.sigma + adj.atom2.sigma) * 0.5
sigma_if_geometric = (adj.atom1.sigma * adj.atom2.sigma) ** 0.5

print(f"{found_14_sigma=}")
print(f"{sigma_if_lorentz=}")
print(f"{sigma_if_geometric=}")
if found_14_sigma != sigma_if_lorentz:
    print("Found 1-4 sigma is not Lorentz-Berthelot!")
if found_14_sigma != sigma_if_geometric:
    print("Found 1-4 sigma is not geometric!")

print("\nApplying hack to rewrite 1-4 interactions with geometric ...")

mod = _patch_14_with_geometric(out)

for adj in mod.adjusts:
    if adj.atom1.name == "C" and adj.atom2.name == "H":
        break

new_14_sigma = adj.type.sigma

if new_14_sigma == sigma_if_lorentz:
    print("Looks like it's still broken \U0001f62C")
elif new_14_sigma == sigma_if_geometric:
    print("Looks fixed! \U0001f60C")
else:
    print("Something is super broken")

Out:

found_14_sigma=2.985
sigma_if_lorentz=2.985
sigma_if_geometric=2.93104077078433
Found 1-4 sigma is not geometric!

Applying hack to rewrite 1-4 interactions with geometric ...
Looks fixed! 😌

Software versions

• Which version of Foyer are you using? (python -c "import foyer; print(foyer.version)") 0.9.1
• Which version of Python (python --version)? 3.8.10
• Which operating system? macOS 11.4

[mattwthompson]

@justinGilmer @daico007 @umesh-timalsina given the severity of this bug (not critical, but also not negligible), could we get a patch release?

[daico007]

@mattwthompson we are planning to do that either today or early tmr

[mattwthompson]

Great - thanks!