ForceField with no Electrostatics tag can be used to make a system

[j-wags]

Describe the bug A ForceField with no Electrostatics tag can be used to create an OpenMM System with arbitrarily-set electrostatics.

To Reproduce

from openforcefield.typing.engines.smirnoff import ForceField
mol=Molecule.from_smiles('O')
ff = ForceField('test_forcefields/tip3p.offxml')
ff.create_openmm_system(mol.to_topology())
from pprint import pprint
pprint(XmlSerializer.serialize(sys))

(?xml version="1.0" ?
 System openmmVersion="7.4.2" type="System" version="1"
   PeriodicBoxVectors
     A x="2" y="0" z="0"/
     B x="0" y="2" z="0"/
     C x="0" y="0" z="2"/
   /PeriodicBoxVectors
   Particles
     Particle mass="15.99943"/
     Particle mass="1.007947"/
     Particle mass="1.007947"/
   /Particles
   Constraints
     Constraint d=".09572000000000001" p1="0" p2="1"/
     Constraint d=".09572000000000001" p1="0" p2="2"/
     Constraint d=".15139006545247014" p1="1" p2="2"/
   /Constraints
   Forces
     **Force alpha="0" cutoff="1" dispersionCorrection="1" 
 ewaldTolerance=".0005" forceGroup="0" ljAlpha="0" ljnx="0" ljny="0" ljnz="0" 
 method="0" nx="0" ny="0" nz="0" recipForceGroup="-1" rfDielectric="78.3" 
 switchingDistance="-1" type="NonbondedForce" useSwitchingFunction="0" 
 version="3"**
       GlobalParameters/
       ParticleOffsets/
       ExceptionOffsets/
       Particles
         Particle eps=".635968" q="-.834" sig=".3150752406575124"/
         Particle eps="0" q=".417" sig="1"/
         Particle eps="0" q=".417" sig="1"/
       /Particles
       Exceptions
         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"/
       /Exceptions
     /Force
   /Forces
 /System)

[mattwthompson]

Could you elaborate on the issue here? It looks like partial charges are populated in the system (unless I'm mis-estimating what that XML means).

Would the solution here just be to error out if trying to call .create_openmm_system with ff['Electrostatics'] missing?

[j-wags]

Ah,sorry to be vague. The issue is with the line:

Force alpha="0" cutoff="1" dispersionCorrection="1" ewaldTolerance=".0005" forceGroup="0" ljAlpha="0" ljnx="0" ljny="0" ljnz="0" method="0" nx="0" ny="0" nz="0" recipForceGroup="-1" rfDielectric="78.3" switchingDistance="-1" type="NonbondedForce" useSwitchingFunction="0" version="3

This has a nonbonded cutoff (1 nm) that was set by the vdWHandler, but is also used for the system's electrostatics. These should, in theory, be separate cutoffs. But the OpenMM system model lumps together vdW and ES into the NonbondedForce, and all we do is have the vdWHandler run first, and then the ElectrostaticsHandler runs second and raises an error if the vdWHandler set the NonBondedForce's cutoff to a different value than it would have used.

Instead, what happens here is that the ElectrostaticsHandler just doesn't exist, so it can't actually verify the cutoff set by the vdWHandler, and the System is made assuming that the electrostatics cutoff is the same as the vdW cutoff.

The desired behavior would be to raise an exception if the ESHandler doesn't exist, but the system has partial charges that will interact.

[mattwthompson]

I ran into this again today (and now better understand the issue) - should test_forcefields/tip3p.offxml be updated to have an <Electrostatics> tag?

[davidlmobley]

I think the answer to your last question is probably yes; and shouldn't we always have an electrostatics tag in the final force field we are applying?

[mattwthompson]

shouldn't we always have an electrostatics tag in the final force field we are applying?

I think @SimonBoothroyd has some niche use cases in which things are split out across separate force fields, but I don't remember too well. I could be confusing it with the need to allow force fields without electrostatics at all.

Fortunately this specific issue is avoided less of an issue with Interchange, since its corresponding electrostatics handler "owns" handlers for each specific charge assignment method and I happened to set the default non-bonded cutoff to 9 Angstroms.

[mattwthompson]

Interchange now attempts to error out in this case.