Missing v-site exceptions with plugins about openff-interchange HOT 5 CLOSED

Your two exception blocks each have 12 exceptions - is it that the number of exceptions are wrong or that the wrong ones are listed? I think the latter, based on what I found below. I need to draw out these exceptions again, I think.

Looks like all the work that went into #906 only hit the combine_nonbonded_forces=True code path; if I try to isolate it to that variable I get what I think is the problematic behavior:

In [1]: import openmm
   ...: from openff.toolkit import Molecule, ForceField, Topology
   ...: from openff.interchange import Interchange
   ...:
   ...: water = Molecule.from_smiles("O")
   ...: water.generate_conformers()
   ...:
   ...: ff = ForceField("tip4p_fb.offxml", load_plugins=True)
   ...: top = Topology.from_molecules([water, water])
   ...: interchange = Interchange.from_smirnoff(force_field=ff, topology=top)
   ...:
   ...: # system = interchange.to_openmm_system(False)
   ...: for combine in [True, False]:
   ...:     system = interchange.to_openmm_system(combine)
   ...:     with open(f"{combine}_water_system.xml", "w") as out:
   ...:         out.write(openmm.XmlSerializer.serialize(system))

In False_water_system.xml, the exceptions are self-consistent but wrong:

	<Forces>
		<Force cutoff="1" energy="4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=(sigma1+sigma2)/2; epsilon=sqrt(epsilon1*epsilon2); " forceGroup="0" method="0" name="vdW force" switchingDistance="-1" type="CustomNonbondedForce" useLongRangeCorrection="0" useSwitchingFunction="0" version="3">
			...
			<Exclusions>
				<Exclusion p1="0" p2="1"/>
				<Exclusion p1="0" p2="2"/>
				<Exclusion p1="1" p2="2"/>
				<Exclusion p1="3" p2="4"/>
				<Exclusion p1="3" p2="5"/>
				<Exclusion p1="4" p2="5"/>
				<Exclusion p1="3" p2="1"/>
				<Exclusion p1="3" p2="2"/>
				<Exclusion p1="7" p2="4"/>
				<Exclusion p1="7" p2="5"/>
				<Exclusion p1="0" p2="3"/>
				<Exclusion p1="3" p2="7"/>
			</Exclusions>
		</Force>
		<Force alpha="0" cutoff="1" dispersionCorrection="1" ewaldTolerance=".0005" exceptionsUsePeriodic="0" forceGroup="0" includeDirectSpace="1" ljAlpha="0" ljnx="0" ljny="0" ljnz="0" method="0" name="Electrostatics force" nx="0" ny="0" nz="0" recipForceGroup="-1" rfDielectric="78.3" switchingDistance="-1" type="NonbondedForce" useSwitchingFunction="0" version="4">
			<Exceptions>
				<Exception eps="0" p1="0" p2="1" q="0" sig="0"/>
				<Exception eps="0" p1="0" p2="2" q="0" sig="0"/>
				<Exception eps="0" p1="1" p2="2" q="0" sig="0"/>
				<Exception eps="0" p1="3" p2="4" q="0" sig="0"/>
				<Exception eps="0" p1="3" p2="5" q="0" sig="0"/>
				<Exception eps="0" p1="4" p2="5" q="0" sig="0"/>
				<Exception eps="0" p1="3" p2="1" q="0" sig="0"/>
				<Exception eps="0" p1="3" p2="2" q="0" sig="0"/>
				<Exception eps="0" p1="7" p2="4" q="0" sig="0"/>
				<Exception eps="0" p1="7" p2="5" q="0" sig="0"/>
				<Exception eps="0" p1="0" p2="3" q="0" sig="0"/>
				<Exception eps="0" p1="3" p2="7" q="0" sig="0"/>
			</Exceptions>
		</Force>

The second molecule's oxygen seems to be the center of these errors; the

The contents of True_water_system.xml seem okay:

				<Exception eps="0" p1="0" p2="1" q="0" sig="1"/>
				<Exception eps="0" p1="0" p2="2" q="0" sig="1"/>
				<Exception eps="0" p1="1" p2="2" q="0" sig="1"/>
				<Exception eps="0" p1="3" p2="4" q="0" sig="1"/>
				<Exception eps="0" p1="3" p2="5" q="0" sig="1"/>
				<Exception eps="0" p1="4" p2="5" q="0" sig="1"/>
				<Exception eps="0" p1="6" p2="1" q="0" sig="1"/>
				<Exception eps="0" p1="6" p2="2" q="0" sig="1"/>
				<Exception eps="0" p1="7" p2="4" q="0" sig="1"/>
				<Exception eps="0" p1="7" p2="5" q="0" sig="1"/>
				<Exception eps="0" p1="0" p2="6" q="0" sig="1"/>
				<Exception eps="0" p1="3" p2="7" q="0" sig="1"/>

(In both cases the atom ordering is O, H, H, O, H, H, V, V)

I bet the problem is that atom indices aren't wired through openff_openmm_particle_map, since it's only the _second_molecule that sees this issue.

I think I need to add tests that at least

• Check that exceptions with virtual sites and without plugins are functionally the same whether or not the non-bonded force is split out
• Do the same with a plugin that modifies the vdW force (in principle this shouldn't change how exceptions are written).
• Do the same with a non-water molecule?

In general,

• It's really important to include more and more systems with multiple molecules with virtual sites in the topology; issues like these don't crop up at all in single-molecule topologies whether or not virtual sites are present
• The logic in creating OpenMM forces (just the non-bonded interactions) is becoming more and more unwieldy over time (up to 1100 lines); I don't think we can afford to slim down any functionality, so it just needs more and more rigorous testing. (Bugs like these, even with the egg on my face, are integral in moving int his direction.)

from openff-interchange.

Ah sorry, to be clear it's that the wrong 12 exceptions are listed. I agree with the testing I was trying to get these energy tests online and also hit this issue, maybe adding something similar where you calculate the energies and forces of a system with stored charges and positions might be useful for this case?

from openff-interchange.

Drawing it out helped - here's two waters, indexed how Interchange organizes particles with OpenMM, and with the exclusions listed out:

This helps make it more clear which exceptions are garbage (3-0, 3-1, 3-2) and by extension which ones are missing. I think there's already a similar test with five-site waters with this sort of analysis, so it might be possible to just run that through with combine_nonbonded_forces=False.

from openff-interchange.

Okay, I found it. The update that ensured all virtual sites were added after all atoms did not cover the split force code path at all. Such a glaring omission can only sneak through via a lack of tests, so that'll be most of the work here.

I even found a simpler way, if accidental, to trigger this behavior. Since the indices are garbled and OpenMM doesn't (by default) add an exception between a pair of atoms where one is already present, just looking at the number of exceptions with multiple waters demonstrates this error. TIP5P, which happened to be the first test I looked at, should produce 10 intramolecular exceptions per water molecule, but a system with 5 waters has slightly fewer than 50 exceptions.

from openff-interchange.

Hopefully fixed with 0.3.23

from openff-interchange.