How to use interchange to save espaloma parameterized system?

[amin-sagar]

Hello. I am trying to save an espaloma parameterized system in gromacs format to run the simulations using gromacs. I posted this issue earlier on espaloma repo https://github.com/choderalab/espaloma/issues/145.

I tried the following script

import os
import torch
import espaloma as esp
from openff.toolkit.topology import Molecule
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout
from openmm.app import Modeller
from simtk import unit
from openmmforcefields.generators import EspalomaTemplateGenerator
from openff.interchange import Interchange
from openff.units.openmm import ensure_quantity

#Load the molecules
molecule = Molecule.from_file("molecule.mol")
## create an Espaloma Graph object to represent the molecule of interest
molecule_graph = esp.Graph(molecule)
# load pretrained model
espaloma_model = esp.get_model("latest")
espaloma_model(molecule_graph.heterograph)
openmm_system = esp.graphs.deploy.openmm_system_from_graph(molecule_graph)
#Load Amber Force Fields 
amber_forcefields = ['amber/protein.ff14SB.xml', 'amber14/tip3pfb.xml']
forcefield = ForceField(*amber_forcefields)
espaloma_generator = EspalomaTemplateGenerator(molecules=molecule, forcefield='espaloma-0.3.1')
openmm_topology = molecule.to_topology().to_openmm()
openmm_positions =  molecule.conformers[0].to_openmm()
#Register Force Field
forcefield.registerTemplateGenerator(espaloma_generator.generator)
#Add water and Ions
modeller = Modeller(openmm_topology, openmm_positions)
modeller.addSolvent(forcefield, model='tip3p', padding=1*nanometer,boxShape='octahedron',ionicStrength=0.15*molar)
system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME,nonbondedCutoff=1.0*nanometer, constraints=HBonds,hydrogenMass=4*amu)

out = Interchange.from_openmm(topology=modeller.topology,system=system, positions=modeller.positions)
# or write to GROMACS files
out.to_gro("out.gro")
out.to_top("out.top")

But I get the following error

 File "/home/amin/mambaforge/envs/espaloma0.3.1/lib/python3.10/site-packages/openff/interchange/components/interchange.py", line 351, in to_gro
    system=_convert(self),
  File "/home/amin/mambaforge/envs/espaloma0.3.1/lib/python3.10/site-packages/openff/interchange/smirnoff/_gromacs.py", line 203, in _convert
    _convert_bonds(molecule, unique_molecule, interchange)
  File "/home/amin/mambaforge/envs/espaloma0.3.1/lib/python3.10/site-packages/openff/interchange/smirnoff/_gromacs.py", line 258, in _convert_bonds
    raise MissingBondError(
openff.interchange.exceptions.MissingBondError: Failed to find parameters for bond with topology indices (2, 3)

I have tried for multiple systems and I always get an error due to missing bond parameters but the simulations run fine. What am I doing wrong? I would be grateful for any suggestions. Amin.

[mattwthompson]

The short answer, and one I expect will get things working for you now, is here. You can still control bond constraints at runtime in the MDP file.

A more thorough solution would involve allowing the GROMACS writer to write topological bonds that are missing physical interactions. The simple cases seem easy enough but I haven't thought through the trickier cases or written up this behavior. I expect that'll happen in the future but I wouldn't count on it happening immediately.

[amin-sagar]

Thanks @mattwthompson I think I understand the origin of the error. However, if I change

system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME,nonbondedCutoff=1.0*nanometer, constraints=HBonds,hydrogenMass=4*amu)

to

system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME,nonbondedCutoff=1.0*nanometer, constraints=None, hydrogenMass=4*amu)

I still get the same error. Does this need to be done with forcefield_kwargs and can't be done by just changing the constraints to None as I tried?

[amin-sagar]

I understand now. The value of rigidWater also needs to be set to False. So, this works.

system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME,nonbondedCutoff=1.0*nanometer,constraints=None,rigidWater=False)

Thanks for the help. Amin.

[mattwthompson]

Hmm ... I'm glad you were able to get it to work but it'd be great if there was a was to use unconstrained bonds and rigid water since TIP3P and its variants are not flexible water models.

[amin-sagar]

That's true. Can someone with more experience in gromacs suggest a way to make waters rigid after topology generation, by adding some restraints?

[mattwthompson]

Usually people use SETTLE or something like it: https://manual.gromacs.org/current/reference-manual/algorithms/constraint-algorithms.html#settle

For example: https://github.com/gromacs/gromacs/blob/e8e7e76946fbe09d952812267580519741796930/share/top/oplsaa.ff/tip3p.itp#L12-L14