vdW handler overwrites v-site charge increments

[SimonBoothroyd]

Describe the bug
The vdW handler currently sets the charges of all atoms to be zero, including those that have had their charge re-distributed by a virtual site.

To Reproduce

from simtk import unit
from simtk.openmm import NonbondedForce

from openff.toolkit.topology import Molecule
from openff.toolkit.typing.engines.smirnoff import ForceField, ParameterList

# Create a FF with a vdW and V-site handler
force_field = ForceField()

constraint_handler = force_field.get_parameter_handler("Constraints")
constraint_handler.add_parameter(
    {"smirks": "[#1:1]-[#8X2H2+0:2]-[#1]", "distance": 0.9 * unit.angstrom}
)
constraint_handler.add_parameter(
    {"smirks": "[#1:1]-[#8X2H2+0]-[#1:2]", "distance": 1.5139 * unit.angstrom}
)

vdw_handler = force_field.get_parameter_handler("vdW")
vdw_handler.add_parameter(
    {
        "smirks": "[*:1]",
        "epsilon": 0.0 * unit.kilojoule_per_mole,
        "sigma": 1.0 * unit.angstrom,
    }
)

virtual_site_handler = force_field.get_parameter_handler("VirtualSites")
virtual_site_handler.add_parameter(
    {
        "smirks": "[#1:1]-[#8X2H2+0:2]-[#1:3]",
        "type": "DivalentLonePair",
        "distance": -0.0106 * unit.nanometers,
        "outOfPlaneAngle": 0.0 * unit.degrees,
        "match": "once",
        "charge_increment1": 1.0552 * 0.5 * unit.elementary_charge,
        "charge_increment2": 0.0 * unit.elementary_charge,
        "charge_increment3": 1.0552 * 0.5 * unit.elementary_charge,
    }
)
virtual_site_handler._parameters = ParameterList(virtual_site_handler._parameters)

# Create an OMM system from the FF
system = force_field.create_openmm_system(Molecule.from_smiles("O").to_topology())

# Print out the charges
force = [
    force for force in system.getForces() if isinstance(force, NonbondedForce)
][0]

for i in range(force.getNumParticles()):
    print(force.getParticleParameters(i)[0])

Output

With the vdW handler added:

0.0 e
0.0 e
0.0 e
-1.0552 e

without the vdW handler added:

0.0 e
0.5276 e
0.5276 e
-1.0552 e

Computing environment (please complete the following information):

• Operating system: OSX
• Output of running conda list: TK 0.9.1

Additional context
Add any other context about the problem here.

[SimonBoothroyd]

cc @j-wags @trevorgokey

[j-wags]

Thanks for reporting this. The seeds of technical debt from the inter-reliance of vdWHandler and ElectrostaticsHandler are coming into bloom 🥀 😞

My instinct is that there should be an error if no ElectrostaticsHandler is present (#716). Does this sound reasonable to you? The correct behavior is achieved if an ElectrostaticsHandler and simple ChargeIncrementHandler are added to the FF:

from simtk import unit
from simtk.openmm import NonbondedForce

from openff.toolkit.topology import Molecule
from openff.toolkit.typing.engines.smirnoff import ForceField, ParameterList

# Create a FF with a vdW and V-site handler
force_field = ForceField()

constraint_handler = force_field.get_parameter_handler("Constraints")
constraint_handler.add_parameter(
    {"smirks": "[#1:1]-[#8X2H2+0:2]-[#1]", "distance": 0.9 * unit.angstrom}
)
constraint_handler.add_parameter(
    {"smirks": "[#1:1]-[#8X2H2+0]-[#1:2]", "distance": 1.5139 * unit.angstrom}
)

vdw_handler = force_field.get_parameter_handler("vdW")
vdw_handler.add_parameter(
    {
        "smirks": "[*:1]",
        "epsilon": 0.0 * unit.kilojoule_per_mole,
        "sigma": 1.0 * unit.angstrom,
    }
)

virtual_site_handler = force_field.get_parameter_handler("VirtualSites")
virtual_site_handler.add_parameter(
    {
        "smirks": "[#1:1]-[#8X2H2+0:2]-[#1:3]",
        "type": "DivalentLonePair",
        "distance": -0.0106 * unit.nanometers,
        "outOfPlaneAngle": 0.0 * unit.degrees,
        "match": "once",
        "charge_increment1": 1.0552 * 0.5 * unit.elementary_charge,
        "charge_increment2": 0.0 * unit.elementary_charge,
        "charge_increment3": 1.0552 * 0.5 * unit.elementary_charge,
    }
)
virtual_site_handler._parameters = ParameterList(virtual_site_handler._parameters)

### NEW CODE HERE ###
ES_handler = force_field.get_parameter_handler("Electrostatics")
charge_increment_handler = force_field.get_parameter_handler("ChargeIncrementModel", {'version':'0.3', 'partial_charge_method': 'formal_charge'})

### END NEW CODE ###

# Create an OMM system from the FF
system = force_field.create_openmm_system(Molecule.from_smiles("O").to_topology())

# Print out the charges
force = [
    force for force in system.getForces() if isinstance(force, NonbondedForce)
][0]

for i in range(force.getNumParticles()):
    print(force.getParticleParameters(i)[0])

yields

0.0 e
0.5276 e
0.5276 e
-1.0552 e

[trevorgokey]

I am not quite following the subtlety in this, but do let me know if it's something caused by vsites, and I'll try to get the fixes in place.

[SimonBoothroyd]

My instinct is that there should be an error if no ElectrostaticsHandler is present (#716). Does this sound reasonable to you? The correct behavior is achieved if an ElectrostaticsHandler and simple ChargeIncrementHandler are added to the FF:

Yeah this was the workaround I ended up finding as well. I think either an error if no handler is found or make the 'vdW' handler an explicit dependency of the v-site handler to be safe.

[mattwthompson]

This appears to be resolved by some combination of the many changes that have happened in the past year. I had to make some modifications to the original reproduction¹² but I believe I preserved the relevant bits:

$ python 885.py
/Users/mattthompson/miniconda3/envs/openff-interchange-env/lib/python3.9/site-packages/openff/toolkit/topology/topology.py:51: TopologyDeprecationWarning: Topology.topology_atoms is deprecated. Use Topology.atoms instead.
  warnings.warn(
0.0 e
0.5276 e
0.5276 e
-1.0552 e

However, this fails when there is no vdW handler present (just deleting the lines involving the vdW handler) - not because Interchange in general requires a vdW handler but because the virtual site code does. There might be a way to implement a workaround but such territory makes me nervous; if particles are going to be added to an OpenMM system, I would much rather their vdW parameters explicitly be specified as zeros via a ParameterHandler with dummy data than implicitly falling back to zeros based on a handler not being found.

from openff.toolkit.topology import Molecule
from openff.toolkit.typing.engines.smirnoff import ForceField, ParameterList
from openff.units import unit
from openmm import NonbondedForce

# Create a FF with a vdW and V-site handler
force_field = ForceField()

# Must create an electrostatics handler for virtual sites to modify charges
electrostatics_handler = force_field.get_parameter_handler("Electrostatics")

# Create a dummy library charge handler, otherwise charges are ambiguous
library_charge_handler = force_field.get_parameter_handler("LibraryCharges")
force_field["LibraryCharges"].add_parameter(
    {
        "smirks": "[*:1]",
        "charge1": "0.0 * elementary_charge",
    }
)

constraint_handler = force_field.get_parameter_handler("Constraints")
constraint_handler.add_parameter(
    {"smirks": "[#1:1]-[#8X2H2+0:2]-[#1]", "distance": 0.9 * unit.angstrom}
)
constraint_handler.add_parameter(
    {"smirks": "[#1:1]-[#8X2H2+0]-[#1:2]", "distance": 1.5139 * unit.angstrom}
)

vdw_handler = force_field.get_parameter_handler("vdW")
vdw_handler.add_parameter(
    {
        "smirks": "[*:1]",
        "epsilon": 0.0 * unit.kilojoule_per_mole,
        "sigma": 1.0 * unit.angstrom,
    }
)

virtual_site_handler = force_field.get_parameter_handler("VirtualSites")
virtual_site_handler.add_parameter(
    {
        "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]",
        "type": "DivalentLonePair",
        "distance": -0.0106 * unit.nanometers,
        "outOfPlaneAngle": 0.0 * unit.degrees,
        "match": "once",
        "charge_increment1": 0.0 * unit.elementary_charge,
        "charge_increment2": 1.0552 * 0.5 * unit.elementary_charge,
        "charge_increment3": 1.0552 * 0.5 * unit.elementary_charge,
    }
)
virtual_site_handler._parameters = ParameterList(virtual_site_handler._parameters)

# Create an OMM system from the FF
system = force_field.create_openmm_system(Molecule.from_smiles("O").to_topology())

# Print out the charges
force = [force for force in system.getForces() if isinstance(force, NonbondedForce)][0]

for i in range(force.getNumParticles()):
    print(force.getParticleParameters(i)[0])

Footnotes

1. I created an ElectrostaticsHandler and also a dummy library charge handler (a ChargeIncrementHandler with partial_charge_method='zeros' would also work) because Interchange enforces that if any ParameterHandler might modify charges an electrostatics handler must be present and there must be some way for atoms to be assigned charges. I believe the existing toolkit does not enforce such requirements and implicitly lets everything remain at zero if unspecified. ↩

2. The original SMIRKS pattern "[#1:1]-[#8X2H2+0:2]-[#1:3]" hit an (explicitly) unsupported case that I didn't take the time to grok. Modifying it to "[#1:2]-[#8X2H2+0:1]-[#1:3]" (and updating the charge increments accordingly) avoided this case, whatever it was. ↩

[mattwthompson]

Double-checked by running the above script with the new builds, and again I get

0.0 e
0.5276 e
0.5276 e
-1.0552 e

which I think it's what's supposed to be reported.