Saving prmtop files from espaloma parameterization

[mmagithub]

Hello,

I am wondering how can we save the espaloma FF parameters (bonded and non-bonded) files, to use in standalone amber simulations, for example, what I need to do after this step:

openmm_system = esp.graphs.deploy.openmm_system_from_graph(molecule_graph,charge_method='am1-bcc')

In order to save ligand.prmtop, ligand.inpcrd, or ligand.off or ligand.frcmod files ?

This is not clear from the documentations.

Thanks,
Marawan

[mikemhenry]

This is a good question, I will have to double-check to see if we support the openff topology object because if we do, we could use https://github.com/openforcefield/openff-interchange to convert to any MD engine that interchange supports

[mmagithub]

this would be awesome, thanks

[mmagithub]

Any update on this request?

[amin-sagar]

I am also trying to do this.
I tried using interchange with 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 an 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 tried using parmed instead using the following script

various imports
import parmed as pmd

#Load the molecules
molecule = Molecule.from_file("Cys-TTZ.mol")
#molecule = pmd.load_rdkit('Cys-TTZ.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)
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()
forcefield.registerTemplateGenerator(espaloma_generator.generator)

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)
system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME,nonbondedCutoff=1.0*nanometer, constraints=None, rigidWater=False)
#out = Interchange.from_openmm(topology=modeller.topology,system=system)
# or write to GROMACS files
#out.to_gro("out.gro")
#out.to_top("out.top")
parmed_structure = pmd.openmm.load_topology(modeller.topology, system=system, xyz=modeller.positions)
#Save Gromacs Files
parmed_structure.save('system.top', overwrite=True)
parmed_structure.save('system.gro', overwrite=True)

#Save amber files
parmed_structure.save('system.parm7', format='amber')
parmed_structure.save('system.rst7', format='rst7')

This allows me to save both gromacs and amber files.
I can run the simulations using the generated files with gromacs, but I am not sure if this is the correct way of doing this.
It would be great to have the opinion of the experts.

Best,
Amin.