openff-toolkit: Incorrect v-site orientation

Describe the bug

The toolkit currently by default sorts the list of parent atom indices when creating a v-site which leads to different orientations for the v-site than were expected.

Presumably this will lead to v-sites being displaced in the wrong direction.

To Reproduce

from openff.toolkit.typing.engines.smirnoff import ForceField
from openmm import unit

force_field = ForceField()

vsite_handler = force_field.get_parameter_handler("VirtualSites")
vsite_handler.add_parameter(
    parameter_kwargs={
        "smirks": "[#6:2]=[#8:1]",
        "name": "EP",
        "type": "BondCharge",
        "distance": 0.7 * unit.nanometers,
        "match": "once",
        "charge_increment1": 0.2 * unit.elementary_charge,
        "charge_increment2": 0.1 * unit.elementary_charge,
        "sigma": 1.0 * unit.angstrom,
        "epsilon": 2.0 / 4.184 * unit.kilocalorie_per_mole,
    }
)

from openff.toolkit.topology import Molecule
molecule = Molecule.from_mapped_smiles("[O:2]=[C:1]=[O:3]")

from openmm import System
omm_system: System
omm_system, topology = force_field.create_openmm_system(
    molecule.to_topology(), return_topology=True
)

for molecule in topology.reference_molecules:
    for v_site in molecule.virtual_sites:
        for particle in v_site.particles:
            print(particle.orientation)

for i in range(omm_system.getNumParticles()):

    if not omm_system.isVirtualSite(i):
        continue

    omm_v_site = omm_system.getVirtualSite(i)

    print(
        tuple(omm_v_site.getParticle(j) for j in range(omm_v_site.getNumParticles()))
    )

Output

Actual

(0, 1)
(0, 2)
(0, 1)
(0, 2)

Expected

(1, 0)
(2, 0)
(1, 0)
(2, 0)

Computing environment (please complete the following information):

  • Operating system
  • Output of running conda list openff-toolkit = 0.10.1

Additional context

About this issue

  • Original URL
  • State: closed
  • Created 3 years ago
  • Reactions: 1
  • Comments: 15 (8 by maintainers)

Commits related to this issue

Most upvoted comments

Yes I agree on all fronts. I was surprised that such a simple example slipped through the existing tests. I think I was hasty in transforming the key since I wanted to keep the “reference” key, but it seems like it was a bad move. It was likely mental carryover from when I was working though all_permutations. I think all of my reservations to the change are only valid for all_permutations, so I think it should be an OK change. Please let me know if something is off after the change and I can help work through any issues 😃

Thanks for the report, Jeff. I checked too and agree it’s broken. I think transforming the key for the match=once case is a mistake.

Yes we likely don’t, but we don’t know if a SMIRKS is completely symmetric so I am not sure how you would detect when to enforce such symmetry constraints