from_smirnoff _without_ charge_from_molecules silently assigns all isormophic molecules the same partial charges

[IAlibay]

Description

User case:

I construct an Interchange object from a Topology with two or more ismorphic Molecules with different partial charges.
In this case the first Molecule has am1bcc charges whilst all the others have gasteiger charges, which differ significantly.

Outcome:

All the Molecules silently get assigned the partial charges from the first Molecule.

Reproduction

from openff.units.openmm import from_openmm
from openmm import NonbondedForce
from openff.toolkit import Molecule, Topology, ForceField
from openff.interchange.components._packmol import solvate_topology_nonwater
from openff.interchange import Interchange
import pytest

solvent = Molecule.from_smiles('c1ccccc1')
solvent.assign_partial_charges(partial_charge_method='gasteiger')
print("solvent charges: ", solvent.partial_charges)

ligand = Molecule.from_smiles('c1ccccc1')
ligand.generate_conformers()
ligand.assign_partial_charges(partial_charge_method='am1bcc')
print("ligand charges: ", ligand.partial_charges)

print("molecules are isomorphic: ", ligand.is_isomorphic_with(solvent))

off_top = Topology.from_molecules(ligand)
solvated_off_top = solvate_topology_nonwater(topology=off_top, solvent=solvent)

ff = ForceField('openff-2.2.0.offxml', 'opc.offxml')
inter = Interchange.from_smirnoff(topology=solvated_off_top, force_field=ff)

for _, c in inter.collections["Electrostatics"].charges.items():
    assert abs(c.m) == pytest.approx(0.13)
    
# Double check this isn't just the collection, and it goes all the way to an openmm System
system = inter.to_openmm_system()
nonbond = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]

for idx in range(system.getNumParticles()):
    c, _, _ = nonbond.getParticleParameters(idx)
    assert abs(from_openmm(c).m) == 0.13

print(all(solvated_off_top.molecule(0).partial_charges == ligand.partial_charges))
print(all(solvated_off_top.molecule(1).partial_charges == ligand.partial_charges))
print(all(solvated_off_top.molecule(1).partial_charges == solvent.partial_charges))

Output

solvent charges:  [-0.062268570782092456 -0.062268570782092456 -0.062268570782092456 -0.062268570782092456 -0.062268570782092456 -0.062268570782092456 0.062268570782092456 0.062268570782092456 0.062268570782092456 0.062268570782092456 0.062268570782092456 0.062268570782092456] elementary_charge
ligand charges:  [-0.13 -0.13 -0.13 -0.13 -0.13 -0.13 0.13 0.13 0.13 0.13 0.13 0.13] elementary_charge
molecules are isomorphic:  True
True
False
True

Software versions

Interchange v0.3.29

[mattwthompson]

For better or worse is the intended behavior - whether or not charge_from_molecules is passed through, the charges on the input topology/molecule objects are completely ignored.

In general this is just how the force field works, partial charges are computed by the definition in the force field and particular molecules can be assigned preset charges by matching the list of molecules in charge_from_molecules. I don't think it's possible to assign isomorphic¹ molecules different charges.

Footnotes

1. Or something like that - I believe the criteria about how similar-looking molecules get the same parameters are baked into Topology.identical_molecule_groups ↩

[IAlibay]

Ah thanks, that makes sense, I totally didn't think through that this was just the partial charges getting reassigned to am1bcc 🤦🏽‍♂️ !

I don't think it's possible to assign isomorphic molecules different charges

That's unfortunate. There are a few cases that come to mind where manually modifying charges for a single molecule could be useful, but I guess it's also something we could do after we convert the Interchange object to the target engine format.

[j-wags]

Yeah, seconding Matt here - this is a known and intentional limitation of our format and ecosystem. We can't assign different charges to isomorphic molecules. If that's the desired behavior then it's up to the user to modify the charges after export to the target engine format.

[davidlmobley]

I agree that this is the intended behavior. Probably the outcome of this is that we should re-check the docs to make sure that this is sufficiently clear in the docs (which I see Matt is working on in #1064 ). However, I'd be opposed to having the toolkit/interchange print something that says it's doing this, because this is the expected behavior.

[IAlibay]

However, I'd be opposed to having the toolkit/interchange print something that says it's doing this, because this is the expected behavior.

@davidlmobley this runs counter to #1048 - from a user perspective, knowing the provenance of how charges were assigned is rather useful (I think I've made my point in that issue but happy to discuss further).