PDBs written without coordinates with OpenEye, sometimes

[lilyminium]

Describe the bug

Sometimes PDB files are written without coordinates with OpenEye. Both times I observed this was from creating a molecule from RDKit.

To Reproduce

(uses MDAnalysis to read and parse the PDB)

from openff.toolkit import Molecule
import MDAnalysis as mda

mol = Molecule.from_smiles("CC")
mol.generate_conformers()
mol.to_file("test.pdb", "PDB")

u = mda.Universe("test.pdb")
u.add_TopologyAttr("elements", u.atoms.types)
rdmol = u.atoms.convert_to("RDKIT")
mol2 = Molecule.from_rdkit(rdmol)
assert len(mol2.conformers)
mol2.to_file("test2.pdb", "PDB")

Output

test.pdb contains coordinates:

> cat test.pdb
ATOM      1  C   UNL     1       0.815  -0.538   0.493  1.00 20.00 
ATOM      2  C   UNL     1       2.155  -0.052  -0.012  1.00 20.00 
ATOM      3  H   UNL     1       0.000   0.001   0.001  1.00 20.00 
ATOM      4  H   UNL     1       0.693  -1.607   0.291  1.00 20.00 
ATOM      5  H   UNL     1       0.730  -0.380   1.572  1.00 20.00 
ATOM      6  H   UNL     1       2.970  -0.591   0.480  1.00 20.00 
ATOM      7  H   UNL     1       2.277   1.017   0.191  1.00 20.00 
ATOM      8  H   UNL     1       2.240  -0.210  -1.091  1.00 20.00 
TER       9      UNL     1 
CONECT    1    2    3    4    5
CONECT    2    6    7    8
MASTER        0    0    0    0    0    0    0    0    8    1    2    0
END

test2.pdb does not:

ATOM      1 C    UNL 
ATOM      2 C    UNL 
ATOM      3 H    UNL 
ATOM      4 H    UNL 
ATOM      5 H    UNL 
ATOM      6 H    UNL 
ATOM      7 H    UNL 
ATOM      8 H    UNL 
TER       9      UNL 
CONECT    1    2    3    4    5
CONECT    2    6    7    8
MASTER        0    0    0    0    0    0    0    0    8    1    2    0
END

If I write it using the RDKitToolkitWrapper though:

ATOM      1  C   UNL     1       0.815  -0.538   0.493  1.00  0.00           C  
ATOM      2  C   UNL     1       2.155  -0.052  -0.012  1.00  0.00           C  
ATOM      3  H   UNL     1       0.000   0.001   0.001  1.00  0.00           H  
ATOM      4  H   UNL     1       0.693  -1.607   0.291  1.00  0.00           H  
ATOM      5  H   UNL     1       0.730  -0.380   1.572  1.00  0.00           H  
ATOM      6  H   UNL     1       2.970  -0.591   0.480  1.00  0.00           H  
ATOM      7  H   UNL     1       2.277   1.017   0.191  1.00  0.00           H  
ATOM      8  H   UNL     1       2.240  -0.210  -1.091  1.00  0.00           H  
CONECT    1    2    3    4    5
CONECT    2    6    7    8
END

Computing environment (please complete the following information):

• Operating system
• Output of running conda list
• Toolkit: 0.16.5 development version (same branch as Allow for lazy loading of other modules with lazy loading #1961)
• OpenEye: 2023.1.0

Additional context

[j-wags]

I've reproduced the issue and it looks like there's something with the residue/metadata info that affects how the PDB is written out. If I clear the hierarchy metadata before writing the PDB the coordinates are written as usual:

from openff.toolkit import Molecule
import MDAnalysis as mda

mol = Molecule.from_smiles("CC")
mol.generate_conformers()
mol.to_file("test.pdb", "PDB")

u = mda.Universe("test.pdb")
u.add_TopologyAttr("elements", u.atoms.types)
rdmol = u.atoms.convert_to("RDKIT")
mol2 = Molecule.from_rdkit(rdmol)
assert len(mol2.conformers)
for atom in mol2.atoms:
    keys = [*atom.metadata.keys()]
    for key in keys:
        del atom.metadata[key]
mol2.to_file("test2.pdb", "PDB")
!cat test2.pdb

returns

ATOM      1 C    UNL     1       0.247  -0.604   0.382  1.00 20.00 
ATOM      2 C    UNL     1      -0.247   0.604  -0.382  1.00 20.00 
ATOM      3 H    UNL     1       1.338  -0.594   0.455  1.00 20.00 
ATOM      4 H    UNL     1      -0.056  -1.527  -0.122  1.00 20.00 
ATOM      5 H    UNL     1      -0.166  -0.613   1.395  1.00 20.00 
ATOM      6 H    UNL     1      -1.338   0.593  -0.455  1.00 20.00 
ATOM      7 H    UNL     1       0.167   0.614  -1.395  1.00 20.00 
ATOM      8 H    UNL     1       0.056   1.527   0.122  1.00 20.00 
TER       9      UNL     1 
CONECT    1    2    3    4    5
CONECT    2    6    7    8
MASTER        0    0    0    0    0    0    0    0    8    1    2    0
END

I'll keep looking into this.

[j-wags]

I can reproduce without MDA - Basically if either insertion_code or chain_id is set on an atom, OETKW doesn't write its coords.

from openff.toolkit import Molecule, OpenEyeToolkitWrapper, RDKitToolkitWrapper

mol_no_mda = Molecule.from_smiles("CC")
mol_no_mda.generate_conformers()
for atom in mol_no_mda.atoms[:3]:
    atom.metadata.update({
        #'residue_name': 'UNL', 
        #'residue_number': 1, 
        'insertion_code': '', 
        #'chain_id': ''
                         })
OpenEyeToolkitWrapper().to_file(mol_no_mda, "test_oe.pdb", file_format='pdb')
!cat test_oe.pdb

yields

ATOM      1  C   UNL     1
ATOM      2  C   UNL     1
ATOM      3  H   UNL     1
ATOM      4  H   UNL     1      -0.056  -1.527  -0.122  1.00 20.00 
ATOM      5  H   UNL     1      -0.166  -0.613   1.395  1.00 20.00 
ATOM      6  H   UNL     1      -1.338   0.593  -0.455  1.00 20.00 
ATOM      7  H   UNL     1       0.167   0.614  -1.395  1.00 20.00 
ATOM      8  H   UNL     1       0.056   1.527   0.122  1.00 20.00 
TER       9      UNL     1 
CONECT    1    2    3    4    5
CONECT    2    6    7    8
MASTER        0    0    0    0    0    0    0    0    8    1    2    0
END

[j-wags]

Can be fixed by replacing insertion_code or chain_id = "" with " ":

from openff.toolkit import Molecule, OpenEyeToolkitWrapper, RDKitToolkitWrapper

mol_no_mda = Molecule.from_smiles("CC")
mol_no_mda.generate_conformers()
for atom in mol_no_mda.atoms[:3]:
    atom.metadata.update({
        #'residue_name': '', 
        #'residue_number': 1, 
        'insertion_code': ' ', 
        #'chain_id': ' '
                         })
OpenEyeToolkitWrapper().to_file(mol_no_mda, "test_oe.pdb", file_format='pdb')
!cat test_oe.pdb

returns

ATOM      1  C   UNL     1       0.247  -0.604   0.382  1.00 20.00 
ATOM      2  C   UNL     1      -0.247   0.604  -0.382  1.00 20.00 
ATOM      3  H   UNL     1       1.338  -0.594   0.455  1.00 20.00 
ATOM      4  H   UNL     1      -0.056  -1.527  -0.122  1.00 20.00 
ATOM      5  H   UNL     1      -0.166  -0.613   1.395  1.00 20.00 
ATOM      6  H   UNL     1      -1.338   0.593  -0.455  1.00 20.00 
ATOM      7  H   UNL     1       0.167   0.614  -1.395  1.00 20.00 
ATOM      8  H   UNL     1       0.056   1.527   0.122  1.00 20.00 
TER       9      UNL     1 
CONECT    1    2    3    4    5
CONECT    2    6    7    8
MASTER        0    0    0    0    0    0    0    0    8    1    2    0
END

[lilyminium]

Thanks for looking into this Jeff, what an odd bug.