I have a molecule with the following SMILES:

sml = "O=C(NCC1=CC(NC(C2=CN=C3C=CC=CN23)=O)=C(C)C=C1)NC4=CC=CC(F)=C4";

The plot looks as follows:

sml // RightComposition[
    Molecule[#, IncludeHydrogens -> True] &,
    MoleculePlot[#, ImageSize -> Large] &

Original molecule showing aromatic bondsIn particular, please note the presence of aromatic bonds.

Now, following the export example at the bottom of the page for PDB files in the help manual, I attempt to export the molecule to a PDB file:

ClearAll[coords, atoms];
{coords, atoms} = sml // RightComposition[
    Molecule[#, IncludeHydrogens -> True] &,
    {MoleculeValue[#, "AtomCoordinates"] &, AtomList /* Map[First]}, Through

    FileNameJoin[{NotebookDirectory[], "mol.pdb"}],
    {{QuantityMagnitude[UnitConvert[coords, "Picometers"]]}, {atoms}},
    {{"AdditionalCoordinates", "AdditionalAtoms"}}

The export happens correctly (or at least so it appears), and I get a PDB file.

But when I load the PDB file, I find that the molecule has lost its aromatic bonds:

"mol.pdb" // RightComposition[
    Import[FileNameJoin[{NotebookDirectory[], #}], {"PDB", "Molecule"}] &,
    MoleculePlot[#, ImageSize -> Large] &

Molecule imported from the PDB fileIn addition, the oxygen atoms appear to have lost their double bonds as well.

What am I doing incorrectly?

Thanks in advance for your help.

  • $\begingroup$ Are you using version 12.1? If so there is a workaround. $\endgroup$
    – Jason B.
    Oct 8, 2020 at 18:45
  • $\begingroup$ @JasonB.yes, I have the latest version, 12.1.1. $\endgroup$ Oct 8, 2020 at 19:01

1 Answer 1

  • PDB files often contain no bond information at all, and when they do they don't contain bond types
  • Perceiving bond types from 3D geometry is a relatively expensive operation, which is fine for small molecules (<100 atoms) but really time-consuming for large proteins
  • The PDB importer for the "Molecule" import element does not do bond-order perception by default, and there isn't a documented option to turn it on.
  • The XYZ importer has a documented option "InferBondTypes" which should be added to the PDB importer, perhaps with an Automatic setting that switches behavior based on the atom count.
  • There is an undocumented option though to turn on the bond order perception for the PDB importer, which can be used in version 12.1
sml = "O=C(NCC1=CC(NC(C2=CN=C3C=CC=CN23)=O)=C(C)C=C1)NC4=CC=CC(F)=C4";
element = {"PDB", "Molecule"};
sml // Molecule // 
    ExportString[#, element]& // 
    ImportString[#, element, "UseOpenBabel" -> True]& //

enter image description here

  • $\begingroup$ Very helpful, thank you! Is there a way to include bond information in the PDB export to being with, so that I won't have to resolve ambiguities later? $\endgroup$ Oct 8, 2020 at 19:37
  • $\begingroup$ I'm not very familiar with the file format to be honest. Reading here it seems that you can specify the presence of bonds but not their type. Using a format like "SDF" makes it easier to round-trip a molecule. $\endgroup$
    – Jason B.
    Oct 8, 2020 at 19:40
  • $\begingroup$ Once again, very helpful to know. I am trying to dock the molecules produced by my WL code using something called AutoDock. Since I am not a chemist, I know almost nothing about these various formats. I will check if AutoDock can accept SDF files. $\endgroup$ Oct 8, 2020 at 19:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.