Openff-toolkit: Unable to export system to CHARMM format

[j-wags]

I'm trying to export a OFF-parameterized ligand to CHARMM format files, and failing. Does anyone have more experience with this?

from openforcefield.topology import Molecule
from openforcefield.typing.engines.smirnoff import ForceField
import parmed

off_mol = Molecule.from_smiles('c1ncccc1OCCCC')
off_mol.generate_conformers()
off_mol.generate_unique_atom_names()
off_topology = off_mol.to_topology()
omm_topology = off_topology.to_openmm(ensure_unique_atom_names=True)
# Load Parsley with H-bond constraints
#off_forcefield = ForceField('openff-1.0.0.offxml')
# Load Parsley without H-bond constraints (assumes CHARMM-GUI will add them later if expected)
off_forcefield = ForceField('openff_unconstrained-1.0.0.offxml')

omm_system = off_forcefield.create_openmm_system(off_topology)
# Convert OpenMM System to a ParmEd structure.
parmed_structure = parmed.openmm.load_topology(omm_topology, omm_system, off_mol.conformers[0])
cps = parmed.charmm.CharmmParameterSet.from_structure(parmed_structure)
parmed.charmm.CharmmParameterSet.write(cps,
                                       top='xxx.rtf',
                                       par='xxx.prm')

---

ParameterError Traceback (most recent call last)
in
----> 1 cps = parmed.charmm.CharmmParameterSet.from_structure(parmed_structure)
2 parmed.charmm.CharmmParameterSet.write(cps,
3 top='p9n.rtf',
4 par='p9n.prm')
5 #parmed_structure.save('p9n.rtf')

~/miniconda3/envs/off-dev/lib/python3.7/site-packages/parmed/charmm/parameters.py in from_structure(cls, struct)
270 """
271 return cls.from_parameterset(
--> 272 ParameterSet.from_structure(struct, allow_unequal_duplicates=False)
273 )
274

~/miniconda3/envs/off-dev/lib/python3.7/site-packages/parmed/parameters.py in from_structure(cls, struct, allow_unequal_duplicates)
296 if not allow_unequal_duplicates and t != dihedral.type:
297 raise ParameterError('Unequal dihedral types defined bewteen '
--> 298 '%s, %s, %s, and %s' % key)
299 break
300 else:

ParameterError: Unequal dihedral types defined bewteen C1, N1, C1, and H1

I've tried manually changing the atomtype names to be distinct in ParmEd:

for idx, atom in enumerate(parmed_structure.atoms):
    atom.atom_type.name = atom.name

But this doesn't work, because ParmEd atom types aren't _strings_, they're _objects_ (so if there are 5 C1 atoms in a molecule, changing one of their types to be called C2 makes them ALL change to a type called C2)

Further debugging:

for omm_atom, pmd_atom in zip(omm_topology.atoms(), parmed_structure.atoms):
    print(omm_atom.id, omm_atom.index, omm_atom.name, pmd_atom.name, pmd_atom.type, pmd_atom.atom_type)

1 0 C1 C1 C1 C1
2 1 N1 N1 N1 N1
3 2 C2 C2 C1 C1
4 3 C3 C3 C1 C1
5 4 C4 C4 C1 C1
6 5 C5 C5 C1 C1
7 6 O1 O1 O1 O1
8 7 C6 C6 C2 C2
9 8 C7 C7 C2 C2
10 9 C8 C8 C2 C2
11 10 C9 C9 C2 C2
12 11 H1 H1 H1 H1
13 12 H2 H2 H1 H1
14 13 H3 H3 H2 H2
15 14 H4 H4 H2 H2
16 15 H5 H5 H3 H3
17 16 H6 H6 H3 H3
18 17 H7 H7 H4 H4
19 18 H8 H8 H4 H4
20 19 H9 H9 H4 H4
21 20 H10 H10 H4 H4
22 21 H11 H11 H4 H4
23 22 H12 H12 H4 H4
24 23 H13 H13 H4 H4

[jchodera]

This is probably a parmed issue to report?

[j-wags]

Probably -- though I'd like to give it a day on this issue tracker in case it's the OpenFF Toolkit's fault, and we can fix it.

[mattwthompson]

I have no experience with CHARMM, but is the issue here that C1, N1, C1, and H1 are being interpreted as "unequal" duplicates? That seems strange

I can get to a typed ethanol structure by a different pathway (with no OpenFF or OpenMM involved) and write it out without triggering this error (although the file written out is garbage, but that's because CHARMM doesn't believe in atom types with names longer than 5 characters).

>>> import mbuild as mb
>>> import foyer
>>> import parmed
>>> ff = foyer.Forcefield(name='oplsaa')
>>> struct = ff.apply(mb.load('CCO', smiles=True).to_parmed())
>>>> # Now we have a parametrized parmed structure
>>> cps = parmed.charmm.CharmmParameterSet.from_structure(struct)     
>>> cps.atom_types
OrderedDict([('OPLS_1', <parmed.topologyobjects.AtomType object at 0x11d409390>)])
>>> parmed.charmm.CharmmParameterSet.write(cps,
...                                        top='xxx.rtf',
...                                        par='xxx.prm')

Okay, so that isn't all that useful. But it at least tells us that the problem is likely related to OpenFF or OpenMM, or the conversion from OpenMM to ParmEd. Here is the source; if I'm reading things correctly the issue is that it already found a dihedral between those four types - not that the printout is just 'C1'- ... but what's actually stored is the pmd.topologyobjects.AtomType _object_, as you pointed out. This seems odd to me, since I would expect there to be duplicate types (more types than unique interactions), not clobbered types (fewer). Maybe something is off with this found_dihed_type_list list in the logic above, I'm not sure.

[j-wags]

Hmm... @mattwthomp could you try the same SMILES as my example? I was
initially using a molecule from file (which caused the error), and ended up
having to make a somewhat elaborate SMILES to replicate the issue.

On Tue, Apr 14, 2020 at 7:46 PM Matt Thompson [e-mail]
wrote:

I have no experience with CHARMM, but is the issue here that C1, N1, C1,
and H1 are being interpreted as "unequal" duplicates? That seems strange

I can get to a typed ethanol structure by a different pathway (with no
OpenFF or OpenMM involved) and write it out without triggering this error
(although the file written out is garbage, but that's because CHARMM
doesn't believe in atom types with names longer than 5 characters).

import mbuild as mb
import foyer
import parmed
ff = foyer.Forcefield(name='oplsaa')
struct = ff.apply(mb.load('CCO', smiles=True).to_parmed())

Now we have a parametrized parmed structure

cps = parmed.charmm.CharmmParameterSet.from_structure(struct)
cps.atom_types
OrderedDict([('OPLS_1', )])
parmed.charmm.CharmmParameterSet.write(cps,
... top='xxx.rtf',
... par='xxx.prm')

Okay, so that isn't all that useful. But it at least tells us that the
problem is likely related to OpenFF or OpenMM, or the conversion from
OpenMM to ParmEd. Here
https://github.com/ParmEd/ParmEd/blob/a8ed12f7d207c81777b29efb7977ef3daa226732/parmed/parameters.py#L297
is the source; if I'm reading things correctly the issue is that it already
found a dihedral between those four types - not that the printout is just 'C1'-
... but what's actually stored is the pmd.topologyobjects.AtomType
object, as you pointed out. This seems odd to me, since I would expect
there to be duplicate types (more types than unique interactions), not
clobbered types (fewer). Maybe something is off with this
found_dihed_type_list list in the logic above, I'm not sure.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
https://github.com/openforcefield/openforcefield/issues/580#issuecomment-613783813,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/AD4STV57BHVAVD2TXVF535TRMUN2DANCNFSM4MIFWHMA
.

[mattwthompson]

c1ncccc1OCCCC writes out files without error using my snippet above. However, I realize now that I was also using a different force field, which adds an unnecessary layer of complexity here

[mattwthompson]

Some updates

• The problem is near the nitrogen and 'c1ncccc1' is sufficient to trigger this issue
• This is the dihedral it's trying to add, but it already finds one with the same periodicity and phase angle but a different k. I had to hunt this down by hand, and can't find the other one, since I don't think there's a clean way to find the source parameters for an interaction in a parametrized system (yet!).

ipdb>  print(t)
<DihedralType; phi_k=0.367, per=2, phase=180.000,  scee=1.000, scnb=1.000>
ipdb>  print(dihedral.type)
<DihedralType; phi_k=2.949, per=2, phase=180.000,  scee=1.000, scnb=1.000>

• The logic here identifies dihedral uniqueness by the tuple of atom type names ((dihedral.atom1.type, dihedral.atom2.type, dihedral.atom3.type, dihedral.atom4.type)), which, despite the name implying they are objects, are actually strings. Elsewhere, there is a mapping between strings and objects.

ipdb>  key
('C1', 'N1', 'C1', 'H1')
ipdb>  (dihedral.atom1.type, dihedral.atom2.type, dihedral.atom3.type, dihedral.atom4.type)
('C1', 'N1', 'C1', 'H1')

In:

unique_atom_types = set([a.atom_type for a in parmed_structure.atoms])
[atom_type.name for atom_type in unique_atom_types]

Out:

['H2', 'H1', 'N1', 'C1']

• Notably, even though @j-wags's debug prints above imply that ParmEd thinks there are several carbon atom types, something happens in the de-duplication process that collapses these down onto a single carbon type.

Could it be that the non-bonded interactions are the same for the different carbons, but dihedrals different on and off the double bond? I would think this could not be the case for an aromatic molecule.

I think this comes down to either an errant dihedral existing when it shouldn't (a problem here) or something unsafe happening in the generation of the atom types (a problem in OpenMM -> ParmEd, I believe).

MWE

```python3
off_mol = Molecule.from_smiles('c1ncccc1')
off_mol.generate_conformers()
off_mol.generate_unique_atom_names()
off_topology = off_mol.to_topology()
omm_topology = off_topology.to_openmm(ensure_unique_atom_names=True)

Load Parsley without H-bond constraints (assumes CHARMM-GUI will add them later if expected)

off_forcefield = ForceField('openff_unconstrained-1.0.0.offxml')

Parametrized and convert OpenMM System to a ParmEd structure

omm_system = off_forcefield.create_openmm_system(off_topology)
parmed_structure = parmed.openmm.load_topology(omm_topology, omm_system, off_mol.conformers[0])

Collect unique parameters

params = parmed.ParameterSet.from_structure(parmed_structure)
print(params.atom_types)

cps = parmed.charmm.CharmmParameterSet.from_structure(parmed_structure) # This will raise an exception
parmed.charmm.CharmmParameterSet.write(cps, top='xxx.rtf', par='xxx.prm')

[j-wags]

Oh, wow. Thanks for diving into this. I hadn't thought that ParmEd might be using nonbonded parameters for "atom type" perception. I wonder if we can make a small, cosmetic change in each nonbonded parameter before exporting to ParmEd to get the atoms "typed" uniquely (if this is the case)

[j-wags]

I think we have this solved in the short term (thanks @mattwthompson!). We can prevent ParmEd from collapsing the "atom types" by changing the nonbonded parameters of each atom by a small amount (1.e-8 nm). MWE [1] is below.

In the long term, we'll want to:

• add a kwarg to ParmEd to prevent it from doing this deduplication
• look at whether the deduplication is based just on nonbonded parameters, or if valence parameters also play a role.
• identify if this sort of thing is a problem outside of just CHARMM conversion, but just a silent error that has been lurking all along

[1]

from openforcefield.topology import Molecule
from openforcefield.typing.engines.smirnoff import ForceField
import parmed
from simtk import openmm, unit

sigma_delta_nm = 1.e-8

off_mol = Molecule.from_smiles('c1ncccc1')
off_mol.generate_conformers()
off_mol.generate_unique_atom_names()
off_topology = off_mol.to_topology()
omm_topology = off_topology.to_openmm(ensure_unique_atom_names=True)

# Load Parsley without H-bond constraints (assumes CHARMM-GUI will add them later if expected)
off_forcefield = ForceField('openff_unconstrained-1.0.0.offxml')


# Parametrized and convert OpenMM System to a ParmEd structure
omm_system = off_forcefield.create_openmm_system(off_topology)
nonbond_force = [force for force in omm_system.getForces() if isinstance(force, openmm.NonbondedForce)][0]

# Add a tiny delta to each radius BEFORE converting to ParmEd
# This makes ParmEd consider these to be unique TYPES when loaded in
for idx in range(nonbond_force.getNumParticles()):
    a, b, c = nonbond_force.getParticleParameters(idx)
    b += sigma_delta_nm * idx * unit.nanometer
    nonbond_force.setParticleParameters(idx, a, b, c)

parmed_structure = parmed.openmm.load_topology(omm_topology, omm_system, off_mol.conformers[0])

# Undo adding the tiny delta (note that this is in Agnstrom, so we multiply delta by 10)
for idx in range(len(parmed_structure.atoms)):
    parmed_structure.atoms[idx].sigma -= sigma_delta_nm * 10 * idx

# Collect unique parameters
params = parmed.ParameterSet.from_structure(parmed_structure)
cps = parmed.charmm.CharmmParameterSet.from_structure(parmed_structure) # This will raise an exception
parmed.charmm.CharmmParameterSet.write(cps, top='xxx.rtf', par='xxx.prm')

[mattwthompson]

The first bullet above is implemented, in part, in https://github.com/ParmEd/ParmEd/pull/1087

[j-wags]

Excellent. I see the above PR was merged. This issue will be closed when that change gets into a ParmEd release. Thanks, Matt!

[jchodera]

@j-wags : To correct this slightly: Since ParmEd releases are decoupled from AmberTools patch-level releases, we need to get that change into an AmberTools patch-level release and then increment the conda-forge version number to automatically build that patch level.

[jaimergp]

This is one thing that I would like to address at some point in the future (having ambertools depend on a parmed package instead of vendoring it, with versions restrained adequately), but we haven't devised a good enough strategy yet!

[j-wags]

Ah, I see. Thanks for the information. What is the normal timescale for ParmEd changes to propagate into AmberTools builds? We have a kludgy workaround for this issue so we don't need this in days, but I also want to make sure we won't be waiting a year for this to get fixed.