Best way to pass espaloma assigned protein-ligand complex system in solution to openff-interchange

[Kamecs]

Dear openforcefield users.

Thank you for your finding of this help.
I want to obtain protein-ligand complex system in solution with espaloma MLFF (https://github.com/choderalab/espaloma).
However, I want to use this system in GROMACS MD, then I want to pass openff.interchange to convert gromacs readable format.

I'am trying to create code to realize it, and following code is currently created:

import sys, os, re
from openff.toolkit.topology import Molecule
from openmm.app import PDBFile, Modeller, ForceField
from openmmforcefields.generators import SystemGenerator
from openff.units.openmm import to_openmm
import openff.units as openff_units
from openmm import unit
from openff.interchange import Interchange

os.environ["INTERCHANGE_EXPERIMENTAL"] = "1"

protein_pdb = PDBFile('protein.pdb')
ligand_molecule = Molecule.from_file('ligand.sdf')
#ligand_molecule.name = "MOL"
#for residue in ligand_molecule.to_topology().to_openmm().residues():
#    residue.name = 'MOL'

modeller = Modeller(protein_pdb.topology, protein_pdb.positions)
modeller.add(
    ligand_molecule.to_topology().to_openmm(),
    to_openmm(ligand_molecule.conformers[0])
)

ff_for_templates = ForceField(
    'amber/ff14SB.xml',
    'amber/tip3p_standard.xml',
    "amber/tip3p_HFE_multivalent.xml"
)

from openmmforcefields.generators import EspalomaTemplateGenerator
espaloma = EspalomaTemplateGenerator(molecules = ligand_molecule)
ff_for_templates.registerTemplateGenerator(espaloma.generator)

modeller.addSolvent(
    ff_for_templates, boxShape = "octahedron",
    model = 'tip3p', padding = 1.0 * unit.nanometers,
    positiveIon = 'Na+', negativeIon = 'Cl-', ionicStrength = 0.15 * unit.molar
)

ic = Interchange.from_openmm(
    topology = modeller.topology,
    positions = modeller.positions,
    box_vectors = modeller.topology.getPeriodicBoxVectors(),
)

I couldn't pass final Interchange.from_openmm part.

...
...
...
RecursionError: maximum recursion depth exceeded while calling a Python object

How can I achieve the desired operation?
Is it really bad to create modeller in the first place?
I attached protein and ligand files.
Archive.zip

Additionally, ligand in my created system is named by UNK.
However, I couln't change residue name of ligand.
Any good way to do that?

I'm sorry to bother you when you're busy.
Best regards.

[mattwthompson]

Unfortunately this is not possible¹. Espaloma, being not really a force field at all, does not organize its interactions in a manner that is easily compatible with GROMACS or Amber which are designed to be used with force fields. If you want to use Espaloma, I recommend sticking with OpenMM unless the original Espaloma developers add new functionality.

Interchange is meant to be less a general molecule mechanics conversion utility and more a means of enabling the use of SMIRNOFF force fields. The existing tools for molecule mechanics conversion (acpype, ParmEd, InterMol, etc.) haven't been used with Espaloma-style parameters and probably won't work.

Footnotes

1. Or, if it is possible, it's not something we've optimized the performance of and won't work for large molecules ↩

[Kamecs]

Dear @mattwthompson

Thank you for your kind and fast reply.

I'm sure that current Espaloma assigned system can not be converted to other file format (e.g. AMBER parm).
Additionally, I also understand that Interchange is created to use SMIRNOFF force fields.

I would like to practice so that I can do things with OpenMM that are comparable to GROMACS + Plumed (Are there any good tutorials out there?).

Best regards.

[mattwthompson]

I'm not a PLUMED user so I don't have any recommendations. There are plenty of resources out there and you're probably fine using whatever is popular.