Exporting all atom hybrid topology to PDB fails

[jchodera]

@dominicrufa : I've been attempting to extract explicit solvent snapshots:

import mdtraj as md
from simtk import unit
import perses
import numpy as np
import os
path = '.'
htf = np.load(os.path.join(path, "htf.npz"), allow_pickle=True)["arr_0"].tolist()
t = md.Trajectory([htf.hybrid_positions / unit.nanometers], htf.hybrid_topology)
t.save('allatom.pdb')

which produces this error

>>> t.save(‘allatom.pdb’)
Traceback (most recent call last):
  File “<stdin>”, line 1, in <module>
  File “/home/chodera/miniconda/envs/perses-new/lib/python3.7/site-packages/mdtraj/core/trajectory.py”, line 1324, in save
    return saver(filename, **kwargs)
  File “/home/chodera/miniconda/envs/perses-new/lib/python3.7/site-packages/mdtraj/core/trajectory.py”, line 1417, in save_pdb
    bfactors=bfactors[i])
  File “/home/chodera/miniconda/envs/perses-new/lib/python3.7/site-packages/mdtraj/formats/pdb/pdbfile.py”, line 672, in __exit__
    self.close()
  File “/home/chodera/miniconda/envs/perses-new/lib/python3.7/site-packages/mdtraj/formats/pdb/pdbfile.py”, line 501, in close
    self._write_footer()
  File “/home/chodera/miniconda/envs/perses-new/lib/python3.7/site-packages/mdtraj/formats/pdb/pdbfile.py”, line 429, in _write_footer
    index1 = atomIndex[atom1]
KeyError: MOL1-C15

Any idea what is going on?

[jchodera]

Attempting to access htf.omm_hybrid_topology is similarly broken:

>>> htf.omm_hybrid_topology
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/chodera/miniconda/envs/perses-new/lib/python3.7/site-packages/perses-0.8.2.dev0-py3.7.egg/perses/annihilation/relative.py", line 2148, in omm_hybrid_topology
    return md.Topology.to_openmm(self._hybrid_topology)
  File "/home/chodera/miniconda/envs/perses-new/lib/python3.7/site-packages/mdtraj/core/topology.py", line 350, in to_openmm
    out.addBond(atom_mapping[a1], atom_mapping[a2], type=bond_mapping[bond.type], order=bond.order)
KeyError: MOL1-C15

[dominicrufa]

Any idea what is going on? KeyError: MOL1-C15

i'm not sure if this is a formatting error, but I'm not sure why the MOL1 prefix is added to the atom name. i presume the error is still thrown if you use t.save_pdb()?

[jchodera]

I don't think it's a formatting error---I think it's an issue with MDTraj's approach to managing bonds, which relies on unique naming conventions for atom names, and I wonder if the issue might be that we have degenerate atom names in our hybrid molecule.

[dominicrufa]

which relies on unique naming conventions for atom names

ah, if that is the case, then the old/new molecules might have unmapped atoms with the same names. even though this is called on each molecule separately, it doesn't know to avoid naming collisions for two distinct molecules.

do you need to save the hybrid as a pdb, or can you save the old/new systems separately? the latter is how we visualize hybrid systems in the visualization PR