Inconsistency between OpenMM simulation run with interchange AMBER and serialized OpenMM system files

[spencercguo]

I am trying to use the OpenFF interchange to create AMBER input files. I am running
something like

topology = solvate_system("1CSA.pdb", ...)
ff = ForceField("modified_forcefield.offxml")
interchange = Interchange.from_smirnoff(force_field=ff, topology=topology)
# save Amber topology and coordinates
interchange.to_prmtop("1CSA.prmtop)
interchange.to_inpcrd("1CSA.inpcrd")
system = interchange.to_openmm_system()
  with open("1CSA.xml', 'w') as output:
      output.write(openmm.XmlSerializer.serialize(system))

and then creating an OpenMM simulation with either

prmtop = AmberPrmtopFile(args.top_file)
inpcrd = AmberInpcrdFile(args.crd_file)
system = prmtop.createSystem(
    nonbondedMethod=PME,
    nonbondedCutoff=1.0 * nanometer,
    constraints=HBonds,)

or

 with open("1CSA.xml", "r") as f:
    system = XmlSerializer.deserialize(f.read())

Then if I create an OpenMM simulation and minimize, I get very different energies (28058442.497487746 kJ/mol
vs -101695.70501193276 kJ/mol). The AMBER system has much higher energy and crashes immediately if I continue the simulation but the deserialized OpenMM simulation works fine.

integrator = LangevinMiddleIntegrator(
    300 * kelvin,
    1 / picosecond,
    2* femtosecond,
)
simulation = Simulation(prmtop.topology, system, integrator)
simulation.context.setPositions(inpcrd.positions)

simulation.minimizeEnergy()

I've attached the original PDB file I used to solvate the system (essentially following this), as well as the resulting prmtop, inpcrd, and openMM XML files. I've also included the forcefield file (a modified version of the OpenFF 2.0) and my conda environment (env.txt).

Archive.zip

[mattwthompson]

Are you able to share a standalone reproduction? I appreciate the work you have put into this but there are unfortunately many details between the PDB file and the prepared system in its various formats. No problem if it's messy or very long! This is complex software and we don't always have the luxury of being succinct.

If not, serializing out to a JSON file is a middle ground.

The first step in debugging this would be figuring out if the problem lies within Interchange (plausible), OpenMM (unlikely), or somewhere in between (my bet). I would like to first run get_summary_data on the prepared system, which requires accessing your solvate_system function.

[spencercguo]

Yep, totally understand! Here are the full solvation and equilibration scripts.

To produce the solvated topologies (prmtop/inpcrd/xml files) I ran the following command
python solvate.py 1CSA.pdb modified_forcefield.offxml 1CSA --solvent water --box-lengths 4 4 4 --method pdbfixer --ionic-strength 0.1
Unfortunately it takes quite a while to run on my machine (I think due to building the interchange needing to recognize all the special parameters in the forcefield file/PDB file).

solvate.py

I also have attached the equilibration script. You can run the script with (the system_file flag is just a workaround to load the serialized Openmm system):
python 1CSA.prmtop 1CSA.inpcrd config.yaml --system-file 1CSA.xml

config.yaml
equil.py

Let me know if you're able to run/reproduce the error (hopefully I'm not missing any files).

[mattwthompson]

Thank you! It might take a day or two for me to get to the bottom of this

[mattwthompson]

I've had a close look at the Amber files you uploaded and something is definitely wrong - hydrogens in water are being given oxygen charges - but I can't isolate it yet. We test similar behavior in similar systems so I would be surprised if something so dire was happening all of the time in Interchange.to_amber, but I can't rule it out.

Could you write out the topology and Interchange objects to disk and upload them here? Something like this around line 178 in solvate.py should do the trick:

    with open("topology.json", "w") as file_topology:
        file_topology.write(topology.to_json())

    with open("interchange.json", "w") as file_interchange:
        file_interchange.write(interchange.model_dump_json())

I'm getting an error trying to reproduce this on my machine:

Traceback (most recent call last):
  File "/Users/mattthompson/Downloads/Archive (1)/solvate.py", line 200, in <module>
    main()
    ~~~~^^
  File "/Users/mattthompson/Downloads/Archive (1)/solvate.py", line 144, in main
    mol = Molecule.from_pdb_and_smiles(args.pdb_file, smi, allow_undefined_stereo=True)
  File "/Users/mattthompson/mamba/envs/toolkit-2136/lib/python3.13/site-packages/openff/toolkit/utils/base_wrapper.py", line 67, in wrapped_function
    value = func(*args, **kwargs)
  File "/Users/mattthompson/mamba/envs/toolkit-2136/lib/python3.13/site-packages/openff/toolkit/topology/molecule.py", line 4792, in from_pdb_and_smiles
    return toolkit.from_pdb_and_smiles(file_path, smiles, allow_undefined_stereo, _cls=cls, name=name)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/mattthompson/mamba/envs/toolkit-2136/lib/python3.13/site-packages/openff/toolkit/utils/rdkit_wrapper.py", line 266, in from_pdb_and_smiles
    raise InvalidConformerError("The PDB and SMILES structures do not match.")
openff.toolkit.utils.exceptions.InvalidConformerError: The PDB and SMILES structures do not match.

using a new environment which is similar to yours

[spencercguo]

Here are the topology and interchange json files. I wonder if it's a problem with PDBfixer, actually? I also tried running with openff-2.3.0-rc3.offxml as the forcefield instead of my custom one and it seems to give the same weird charges for water in the prmtop.

interchange.zip

[mattwthompson]

Thanks for the files! I can load them up and am starting to look through them.

I'm also pulling one of our resident PDB experts @Yoshanuikabundi into the fold here - could you describe how you prepared this file? Cyclic proteins are non-trivial to load and there may be better ways to load this one into OpenFF tooling.

[spencercguo]

The PDB file is of cyclosporin from the PDB (solution NMR) and it's been lightly minimized/prepped in Maestro, I believe. I was having trouble figuring out the best way of loading/handling non canonical amino acids, which is how I ended up with the strange method that is in the script.

[Yoshanuikabundi]

We're currently polishing up a new way of loading PDB files that's designed to be easy to use with NCAAs. This PDB file has a couple of features that make it a bit difficult, including that Maestro has... rearranged... its residue order, but it has enough CONECT records that we should be able to handle it. I'll fold it into the test suite, fix the issues, and get back to you!

[spencercguo]

We're currently polishing up a new way of loading PDB files that's designed to be easy to use with NCAAs. This PDB file has a couple of features that make it a bit difficult, including that Maestro has... rearranged... its residue order, but it has enough CONECT records that we should be able to handle it. I'll fold it into the test suite, fix the issues, and get back to you!

Awesome, I had just seen that prototype you linked, but haven't had the chance to try it out. Definitely seems like it'll make life easier handling NCAAs.

[mattwthompson]

In an effort to keep this a little more organized, I'm separating this out into two parts:

1. PDB loading, which Josh is working on in some capacity (I won't speak for them on timelines, etc.) and potential errors in the file itself
2. Amber vs. OpenMM energies: this is clearly a bug and I'm working on it in a separate child issue. Please follow there for updates on that part!

[mattwthompson]

@spencercguo please upgrade to Interchange 0.4.10 which will fix the Amber issue

[spencercguo]

@spencercguo please upgrade to Interchange 0.4.10 which will fix the Amber issue

Awesome, thanks for helping check this out and getting the fix!