Combining SMIRNOFF interchanges with isomorphic molecules with different charges uses only the charges from the latter molecule

[Yoshanuikabundi]

Description

I was looking to see if @IAlibay's desired workflow in #1059 could be accomplished by creating two separate interchanges and then combining them to get a single interchange with different charges for isomorphic molecules, and it does not. I think this is surprising - in my head, an interchange associates each atom in the topology with its own parameters, and the fact that parameters are deduplicated is like an optimization. I've already unambiguously identified which atoms should have which charges, and so to lose that information for that optimization is surprising.

Reproduction

Modified version of Irfan's code:

from openff.toolkit import Molecule, Topology, ForceField
from openff.interchange.components._packmol import solvate_topology_nonwater
from openff.interchange import Interchange
from openff.units import unit
import numpy as np

solvent = Molecule.from_smiles('c1ccccc1')
solvent.assign_partial_charges(partial_charge_method='gasteiger')

ligand = Molecule.from_smiles('c1ccccc1')
ligand.generate_conformers()
ligand.assign_partial_charges(partial_charge_method='am1bcc')

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)
solvent_off_top = Topology.from_molecules([solvated_off_top.molecule(i) for i in range(off_top.n_molecules, solvated_off_top.n_molecules)])

ff = ForceField('openff-2.2.0.offxml', 'opc.offxml')
ligand_interchange = Interchange.from_smirnoff(topology=off_top, force_field=ff, charge_from_molecules=[ligand])
inter_lig_charges = [c.m_as(unit.elementary_charge) for c in ligand_interchange.collections["Electrostatics"].charges.values()] * unit.elementary_charge

solvent_interchange = Interchange.from_smirnoff(topology=off_top, force_field=ff, charge_from_molecules=[solvent])
inter_sol_charges = [c.m_as(unit.elementary_charge) for c in solvent_interchange.collections["Electrostatics"].charges.values()] * unit.elementary_charge

inter = ligand_interchange.combine(solvent_interchange)
inter_charges = [c.m_as(unit.elementary_charge) for c in inter.collections["Electrostatics"].charges.values()] * unit.elementary_charge

print("ligand charges in ligand interchange: ", inter_lig_charges)
print("charges in solvent interchange: ", inter_sol_charges)
print("charges in combined interchange: ", inter_charges)

assert not np.isclose(inter_charges, np.append(inter_lig_charges, inter_sol_charges)).all()

Output

molecules are isomorphic:  True
ligand charges in ligand interchange:  [-0.1301600026587645 -0.1300999956826369 -0.1300999956826369 -0.1300999956826369 -0.1300999956826369 -0.1300999956826369 0.13010999684532484 0.13010999684532484 0.13010999684532484 0.13010999684532484 0.13010999684532484 0.13010999684532484] elementary_charge
charges in solvent interchange:  [-0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489] elementary_charge
charges in combined interchange:  [-0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 -0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489 0.06175815314054489] elementary_charge

Software versions

Interchange v0.4.0beta2

[mattwthompson]

I'm still iffy on whether or not this is a valid use case, but I agree that the behavior you're seeing is surprising

[j-wags]

I (7/10) think this behavior is valid - Generally, it should be fine to combine interchanges created by two different SMIRNOFF FFs, without either one of them having its assigned parameters overwritten.

[mattwthompson]

I should expand a little bit on my views -

• Using different charge models for (approximately) identical molecules: walking a tightrope between funky and invalid
• Getting different charges on isomorphic molecules from one force field: not valid
• Combining the results of (intentionally) applying different charge models to topologies which happen to contain an isomorphic molecule(s): quite funky, but not wholly invalid - and the charges on the output should reflect the charges of each input

[Yoshanuikabundi]

I could imagine someone running a simulation with a deliberately conformation-dependent charge model on two different conformations of the same molecule, on a timescale where those conformations don't interconvert. Funky yes, but pretty reasonable thing to do, definitely valid. So I don't think it's completely obvious that different charges on isomorphic molecules from one force field is invalid? I agree it's currently impossible to write such a force field in SMIRNOFF without plugins.

I could also imagine different charges for the same molecule if those molecules are isolated in different environments for the entire simulation (gas and liquid phases separated by some membrane?). Or maybe alchemic mutation on a subset of copies of a molecule. In Martini 2 you need 10% of your water beads to have different LJ parameters to the rest so that it doesn't freeze at like 20 Celsius, I know coarse grained models are out of scope for now but my point is that people do all sorts of weird things to make their force fields work.

I can also imagine maybe one day this helps someone debug a problem with a different combined system, or make a theoretically indefensible but numerically correct simulation work or something weird like that. As a general rule, in a perfect world I think combined interchanges should have the same parameters for the same atoms as the original interchanges, and that users shouldn't have to think about any edge cases for this conceptually pretty well defined operation. Not making any claims about whether that makes the effort to fix this bug worth it.

[mattwthompson]

So I don't think it's completely obvious that different charges on isomorphic molecules from one force field is invalid?

This case is invalid if and only if it comes from a single force field, and I think that's the key difference in perspectives here; lots of stuff is possible, but fewer things are really well-defined in force fields

For excessively exotic use cases, the user can also just modify charges after the fact (needs to be documented better sometime #1071)

[Yoshanuikabundi]

When you say "single force field", do you mean, like, a single conceptual force field, or a single SMIRNOFF force field? For example I could imagine having a liquid phase force field and a gas phase force field that are each a single SMIRNOFF force field, but I think conceptually of the combination of the two as a single force field and sometimes use them together in the same simulation.

[mattwthompson]

I mean a single object

[mattwthompson]

Slightly slimming down the first code snippet (okay, well, ... there's only so much that can be stripped out):

In [1]: from openff.toolkit import Molecule, Topology, ForceField, Quantity
   ...: from openff.interchange import Interchange
   ...: from openff.units import unit
   ...:
   ...: solvent = Molecule.from_smiles('c1ccccc1')
   ...: solvent.assign_partial_charges(partial_charge_method='gasteiger')
   ...:
   ...: ligand = Molecule.from_smiles('c1ccccc1')
   ...: ligand.assign_partial_charges(partial_charge_method='am1bcc')
DEPRECATION: --no-python-version-warning is deprecated. pip 25.1 will enforce this behaviour change. A possible replacement is to remove the flag as it's a no-op. Discussion can be found at https://github.com/pypa/pip/issues/13154

In [2]: sage = ForceField("openff-2.2.0.offxml")
/Users/mattthompson/micromamba/envs/interchange-examples-env/lib/python3.12/site-packages/smirnoff99frosst/smirnoff99frosst.py:11: UserWarning: Module openff was already imported from None, but /Users/mattthompson/software/openff-interchange is being added to sys.path
  from pkg_resources import resource_filename

In [3]: int1 = sage.create_interchange(solvent.to_topology(), charge_from_molecules=[solvent])

In [4]: int2 = sage.create_interchange(solvent.to_topology(), charge_from_molecules=[ligand])

In [5]: int1['Electrostatics'].charges
Out[5]:
{TopologyKey with atom indices (0,): -0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (1,): -0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (2,): -0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (3,): -0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (4,): -0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (5,): -0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (6,): 0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (7,): 0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (8,): 0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (9,): 0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (10,): 0.06175815314054489 <Unit('elementary_charge')>,
 TopologyKey with atom indices (11,): 0.06175815314054489 <Unit('elementary_charge')>}

In [6]: int2['Electrostatics'].charges
Out[6]:
{TopologyKey with atom indices (0,): -0.1301600026587645 <Unit('elementary_charge')>,
 TopologyKey with atom indices (1,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (2,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (3,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (4,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (5,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (6,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (7,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (8,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (9,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (10,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (11,): 0.13010999684532484 <Unit('elementary_charge')>}

In [7]: int1.combine(int2)['Electrostatics'].charges
/Users/mattthompson/software/openff-interchange/openff/interchange/operations/_combine.py:52: InterchangeCombinationWarning: Interchange object combination is complex and likely to produce strange results outside of a limited set of use cases it has been tested in. Any workflow using this method is not guaranteed to be suitable for production or stable between versions. Use with extreme caution and thoroughly validate results!
  warnings.warn(
/Users/mattthompson/software/openff-interchange/openff/interchange/operations/_combine.py:124: InterchangeCombinationWarning: Setting positions to None because one or both objects added together were missing positions.
  warnings.warn(
Out[7]:
{TopologyKey with atom indices (0,): -0.1301600026587645 <Unit('elementary_charge')>,
 TopologyKey with atom indices (1,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (2,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (3,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (4,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (5,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (6,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (7,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (8,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (9,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (10,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (11,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (12,): -0.1301600026587645 <Unit('elementary_charge')>,
 TopologyKey with atom indices (13,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (14,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (15,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (16,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (17,): -0.1300999956826369 <Unit('elementary_charge')>,
 TopologyKey with atom indices (18,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (19,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (20,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (21,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (22,): 0.13010999684532484 <Unit('elementary_charge')>,
 TopologyKey with atom indices (23,): 0.13010999684532484 <Unit('elementary_charge')>}