openforcefield/openff-forcefields

Detected Incomplete Exceptions and Unassigned Bond Parameters when converting to GROMACS system using ParmEd

ajfriedman22 opened this issue · 1 comments

A warning was generated when attempting to convert the generated OpenMM system to GROMACS. Additionally, several bonds do not have parameter terms assigned (see AD_OpenFF.top, line 84 for atoms 2-7, 94 for atoms 5-13, etc). The code used is below as well as the generated error and the PDB file. I have also included the itp file generated using CGennFF for the same molecule and the topology file which was generated with ParmEd.
AD_OpenFF.zip

#Loading OpenFF Forcefield and parmed
from openff.toolkit.typing.engines.smirnoff.forcefield import ForceField
from simtk.openmm.app import PDBFile
import parmed
from simtk import unit
import numpy as np

parsley = ForceField('openff_unconstrained-1.3.0.offxml')

#Loading up an OpenFF Topology from amorphadiene's SMILES
from openff.toolkit.topology import Molecule, Topology

AD = Molecule.from_smiles('CC1CCC2C(C)CCC(C2C1)C(C)=C', allow_undefined_stereo=True)

#Import PDB file
pdbfile = PDBFile('AD_processed.pdb')

#Generate topology
top_ff = Topology.from_molecules(1 * [AD])
top_mm = pdbfile.topology

#Set box vectors for OpenFF topology
top_ff.box_vectors = np.array([4, 4, 4]) * unit.nanometer

#Create an OpenMM System
openmm_sys = parsley.create_openmm_system(top_ff)

#Convert OpenMM System to a ParmEd structure.
parmed_structure = parmed.openmm.load_topology(top_mm, openmm_sys, pdbfile.positions)

#Export GROMACS files.
parmed_structure.save('AD_OpenFF.top', overwrite=True)
parmed_structure.save('AD_OpenFF.gro', overwrite=True)
/home/afriedman/doc/miniconda3/lib/python3.8/site-packages/parmed/openmm/topsystem.py:451: OpenMMWarning: Detected incomplete exceptions. Not supported.
  warnings.warn('Detected incomplete exceptions. Not supported.', OpenMMWarning)

I think the problem here arises from how the different OpenMM Topology objects were created; openmm_sys is created from top_ff, not pdbfile.topology, and something is causing them to become out of sync. Even though they're the same type and they're similar enough that you can use the positions of one to set the other, something about them is different enough to cause ParmEd to get confused. I don't know why that's the case, but that seems to be the root of the problem. I can see two ways to fix this, below. In both cases, I recommend using top_ff.to_openmm() to make the OpenMM topology passed to ParmEd, since this is the same object that was passed through the system generation step. Locally, these write out GROMACS files that may or may not be accurate, but pass the eye check (no missing bond parameters). I don't think there's a way for us to enforce that users use compatible objects in the create_openmm_system and load_topology calls since ParmEd is a downstream package and necessary for this use case (for now).

1. Don't use the PDB file at all

This has the benefit of using only one source of data, not relying on the PDB file at all. However, this would make it harder to run multi-molecule systems, since positions of one conformer of one molecule are hard to expand out to other systems.

mol = Molecule.from_smiles('CC1CCC2C(C)CCC(C2C1)C(C)=C', allow_undefined_stereo=True)
mol.generate_conformers(n_conformers=1)
positions = mol.conformers[0].in_units_of(unit.nanometer)

top = Topology.from_molecules(1 * [mol])

...
openmm_sys = parsley.create_openmm_system(top)

parmed_structure = parmed.openmm.load_topology(top.to_openmm(), openmm_sys, positions)
...

2. Construct the OpenFF topology from the OpenMM topology

This is most similar to the GROMACS example. It has the benefit of working with multiple molecules, i.e. if you can make a PDB file of 100 of these molecules with something like PACKMOL, the new position data is there.

mol = Molecule.from_smiles('CC1CCC2C(C)CCC(C2C1)C(C)=C', allow_undefined_stereo=True)

pdbfile = PDBFile('AD_processed.pdb')
top = Topology.from_openmm(pdbfile.topology, unique_molecules=[mol])

openmm_sys = parsley.create_openmm_system(top)

parmed_structure = parmed.openmm.load_topology(top.to_openmm(), openmm_sys, pdbfile.positions)

I'm going to close this issue since it doesn't seem to be a problem of missing parameters, but feel free to ask follow-ups if something else comes up.