Incorrect bond perception from 3D

[ldgibson]

• I believe this to be a bug with Open Babel
• This is a feature request

Environment Information

Open Babel version: 2.4.1 (installed using conda install -c openbabel openbabel)
Operating system and version: MacOS 10.12.6

Expected Behavior

Expected SMILES string: [Li]OC(=O)OCC[O+]1CCOC1=O

Actual Behavior

Generated SMILES string: [Li]OC(=O)OC[CH2].C1(=O)OCCO1

Steps to Reproduce

After copying the structure information into structure.xyz, the following code should reproduce the bug.

import pybel

mol = next(pybel.readfile('xyz', 'structure.xyz'))
smiles = mol.write('smiles').split('\t')[0]
print(smiles)

XYZ file contents:

21
    Energy:      46.5530552
Li        -1.05918       -6.48375        3.03497
C         -3.67179       -3.25569        4.73488
O         -2.88169       -4.12988        5.05149
O         -3.95901       -6.03517        3.39980
O         -4.97821       -3.46526        4.48731
C         -5.26948       -5.59870        3.34006
C         -5.65864       -4.68124        4.51417
H         -5.94439       -6.48245        3.31165
H         -5.41690       -5.03354        2.39534
H         -6.75045       -4.49118        4.40348
H         -5.51590       -5.18271        5.49492
C         -3.47339       -6.91364        4.31571
O         -2.16776       -7.23753        4.29134
O         -3.32406       -1.96706        4.65992
O         -4.20346       -7.42376        5.15235
C         -4.43330       -1.22706        4.32170
C         -5.58838       -2.24603        4.27319
H         -4.26286       -0.76996        3.32316
H         -4.62671       -0.43724        5.07871
H         -6.09654       -2.21804        3.28414
H         -6.32114       -2.05494        5.08845

[ghutchis]

I see your point, but I haven't had many people clamoring to see trivalent O+ species.

The problem with the bond perception code is that it has to serve many masters - the most frequent being people reading PDB or CIF files, where there's atomic disorder.

So right now, the code restricts some elements to a maximum valence (2 for oxygen) and drop the longest bond - here the C-O.

My guess is that increasing valence to 3 for oxygen is going to have far more problems outside your use case. Your profile mentions MD - my strong suggestion is to either write your own bond perception code (for single bonds it's pretty easy) or track the bond topology through the MD.

[ldgibson]

I appreciate your fast reply.

As far as writing my own bond perception code, I track bonds using distance cutoffs during simulation and from there I build contact matrices. So I essentially know what the topology of each molecule is, but I still preferred to use SMILES strings as there are many tools for visualization. Do you know if it is possible for me to instead somehow pass my system's contact matrix to openbabel and have it generate the SMILES string? Or am I simply better off building my own visualization tools for my tracked contact matrices?

Thank you for the help!

[ghutchis]

No, no need to build your own tools. That's the whole point of Open Babel and Avogadro.

There are a few options, but the easiest is probably to read in the XYZ file, delete all the bonds, and then create them as needed with mol.OBMol.AddBond(beginIndex, endIndex, bondOrder)

Keep in mind that to delete the bonds, you'll want to iterate through them, add to a list and then delete them through the list:

    bonds = []
    for bond in ob.OBMolBondIter(mol.OBMol):
        bonds.append(bond)  # if you delete the bond, the OBMolBondIter will get messed up
    for bond in bonds:
        mol.OBMol.DeleteBond(bond)