Virtual site parameters can clash if looked up using only SMIRKS

[lilyminium]

Description

I was testing out some match="once" combinations with TrivalentLonePair and ran into an unexpected key/potential mapping in Interchange.

Reproduction

Please include a minimally reproducing example of this bug.

from openff.toolkit import Molecule, ForceField
import numpy as np
from openff.units import unit

ff = ForceField("openff-2.1.0.offxml")
handler = ff.get_parameter_handler("VirtualSites")
handler.add_parameter({
  "type": "TrivalentLonePair",
  "smirks": "[#7:1](-[*:2])(-[*:3])-[*:4]",
  "name": "EP1",
  "distance": 0.5 * unit.angstrom,
  "charge_increment": np.array([1, 0, 0, 0]) * unit.elementary_charge,
  "sigma": 0.0 * unit.angstrom,
  "epsilon": 0.0 * unit.kilojoules_per_mole,
  "match": "once"
})
handler.add_parameter({
  "type": "TrivalentLonePair",
  "smirks": "[#7:1](-[*:2])(-[*:3])-[*:4]",
  "name": "EP2",
  "distance": -1 * unit.angstrom,
  "charge_increment": np.array([5, 0, 0, 0]) * unit.elementary_charge,
  "sigma": 0.0 * unit.angstrom,
  "epsilon": 0.0 * unit.kilojoules_per_mole,
  "match": "once"
})

mol = Molecule.from_smiles("N")
ic = ff.create_interchange(mol.to_topology())

Output

# print(ic.collections["VirtualSites"].potentials)
{PotentialKey associated with handler 'VirtualSites' with id '[#7:1](-[*:2])(-[*:3])-[*:4] EP1 once': Potential(parameters={'distance': <Quantity(0.5, 'angstrom')>, 'outOfPlaneAngle': None, 'inPlaneAngle': None}, map_key=None),
 PotentialKey associated with handler 'VirtualSites' with id '[#7:1](-[*:2])(-[*:3])-[*:4] EP2 once': Potential(parameters={'distance': <Quantity(0.5, 'angstrom')>, 'outOfPlaneAngle': None, 'inPlaneAngle': None}, map_key=None)}

# v1, v2 = ic.collections["VirtualSites"].potentials.values()
# print(v1 is v2)
False
# print(v1 == v2)
True

# ic.collections["Electrostatics"].charges
{TopologyKey with atom indices (0,): 0.9872900098562241 <Unit('elementary_charge')>,
 TopologyKey with atom indices (1,): 0.337569996714592 <Unit('elementary_charge')>,
 TopologyKey with atom indices (2,): 0.337569996714592 <Unit('elementary_charge')>,
 TopologyKey with atom indices (3,): 0.337569996714592 <Unit('elementary_charge')>,
 VirtualSiteKey with atom indices None: -1.0 <Unit('elementary_charge')>,
 VirtualSiteKey with atom indices None: -1.0 <Unit('elementary_charge')>}

Despite Interchange picking up both EP1 and EP2, EP2 appears to get assigned the same parameters as EP1. If I reverse the order of EP2 and EP1 in defining the force field, the reverse happens -- EP1 gets assigned the same parameters as EP2.

Software versions

• Which operating system and version are you using?
• Mac, main
• How did you install Interchange? pip install -e .
• What is the output of running conda list?

$ micromamba list openff
 Name                   Version    Build         Channel
───────────────────────────────────────────────────────────────
  openff-amber-ff-ports  0.0.4      pyhca7485f_0  conda-forge
  openff-forcefields     2024.03.0  pyhca7485f_0  conda-forge
  openff-models          0.0.5      pyh1a96a4e_0  conda-forge
  openff-nagl            0.2.2      pyhd8ed1ab_0  conda-forge
  openff-nagl-models     0.0.2      pyhd8ed1ab_0  conda-forge
  openff-qcsubmit        0.50.2     pyhd8ed1ab_0  conda-forge
  openff-toolkit         0.15.2     pyhd8ed1ab_0  conda-forge
  openff-toolkit-base    0.15.2     pyhd8ed1ab_0  conda-forge
  openff-units           0.2.1      pyh1a96a4e_0  conda-forge
  openff-utilities       0.1.8      pyh1a96a4e_0  conda-forge

[mattwthompson]

Okay, so you have an ammonia and you're trying to apply two different virtual sites using TrivalentLonePair. They have the same SMIRKS pattern but different names, geometries, and charges (neither have vdW parameters). AM1BCC gives ammona roughly -1e on the nitrogen and +1/3e on the hydrogens, so with these values you'd expect to end up with a +5e on the nitrogen and two virtual sites with -1e and -5e. Similarly, you'd expect to have virtual sites named EP1 and EP2 with positive and negative distances, respectively. But what you get is separate names with all of the parameters of EP2 mushed onto EP1, like EP2 was never seen.

In [89]: [key.name for key in ic["VirtualSites"].key_map]
Out[89]: ['EP1', 'EP2']

In [90]: [pot.parameters["distance"] for pot in ic["VirtualSites"].potentials.values()]
Out[90]: [<Quantity(0.5, 'angstrom')>, <Quantity(0.5, 'angstrom')>]

In [91]: pprint([*ic["Electrostatics"].potentials.values()])
[Potential(parameters={'charge': <Quantity(-1.01270999, 'elementary_charge')>}, map_key=None),
 Potential(parameters={'charge': <Quantity(0.337569997, 'elementary_charge')>}, map_key=None),
 Potential(parameters={'charge': <Quantity(0.337569997, 'elementary_charge')>}, map_key=None),
 Potential(parameters={'charge': <Quantity(0.337569997, 'elementary_charge')>}, map_key=None),
 Potential(parameters={'charge_increments': <Quantity([1 0 0 0], 'elementary_charge')>}, map_key=None),
 Potential(parameters={'charge_increments': <Quantity([1 0 0 0], 'elementary_charge')>}, map_key=None)]

In [92]: [charge.m for charge in ic["Electrostatics"].charges.values()]
Out[92]: [0.9872900098562241, 0.337569996714592, 0.337569996714592, 0.337569996714592, -1.0, -1.0]

This is definitely wrong. That the keys aren't wrong implies to me it's something going wrong in the logic that decides if there should be a new parameter or not, which shouldn't happen since it gets both results from SMARTS-matching.

Incidentally, I found that the position of each is getting assigned incorrectly. OpenMM should be smart enough to fix this if the parameters are correct, so I can often wave my hands and claim it's just a visualization quirk. But in this case the parameters themselves appear to also get the value of EP1:

			<LocalCoordinatesSite p1="0" p2="1" p3="2" p4="3" pos1="-.049999999999999996" pos2="0" pos3="0" wo1="1" wo2="0" wo3="0" wo4="0" wx1="-1" wx2=".3333333333333333" wx3=".3333333333333333" wx4=".3333333333333333" wy1="-1" wy2="1" wy3="0" wy4="0"/>
		</Particle>
		<Particle mass="0">
			<LocalCoordinatesSite p1="0" p2="1" p3="2" p4="3" pos1="-.049999999999999996" pos2="0" pos3="0" wo1="1" wo2="0" wo3="0" wo4="0" wx1="-1" wx2=".3333333333333333" wx3=".3333333333333333" wx4=".3333333333333333" wy1="-1" wy2="1" wy3="0" wy4="0"/>

Might be related, might not. With any luck, the charge and distance parameters are mis-assigned from only the virtual site code and the actual logic of getting charges and positions are not incorrect.

[mattwthompson]

This is probably where things are getting funky:

https://github.com/openforcefield/openff-interchange/blob/v0.3.24/openff/interchange/smirnoff/_virtual_sites.py#L145

This has caused me several headaches in the past and I've advocated for it to be changed (openforcefield/openff-toolkit#1363) but I think with virtual sites it's valid for two fundamentally different parameters to have the same SMIRKS patterns.

The lookup magic is poorly-suited for this case, so I think I'll have to temporarily carry along the entire VirtualSiteType.

[mattwthompson]

Oh also @lilyminium I'm curious if this test cases was really a MWE or something closer to scientific interest. Covering this with a unit test(s) should fix the bug, but working it into an example increases the likelihood it really works (i.e. can run a few timesteps in addition to having the appearance of being parameterized correctly).

In other words, is there a practical use case of these sorts of virtual sites, being two particles at difference distances and/or charges on the same orientation atom(s)?

[lilyminium]

This particular MWE wasn't of scientific interest, but I think I would have liked this behaviour for easier multi-vsites on the C-Br bond (which I currently implemented using both C and Br as parent atoms, but it would be easier to use Br as the parent for both) -- I was testing things out in relation to openforcefield/standards#64