How to optimize the molecular geometry under PyMol?
Disclaimer: I don't typically use PyMol, so I only made a quick look at the Optimize.py script.

You're right that the problem is with the bond orders. The Optimize.py script uses a PDB file as an intermediate representation that it sends to Open Babel for force field optimization and other properties.

The problem is that PDB does not support bond orders. So yes, you don't get double-bond information back, and I can't verify if Open Babel is getting bond orders from PyMol and Optimize.py.

I looked through the PyMol commands and it seems as if there's a way to write an MDL Molfile (which would have bond orders) but not read one - PDB is the only option for reading and writing in a PyMol script.

My suggestions:


Either push the PyMol developers to support other formats (e.g., Sybyl .mol2 or MDL Mol) in the scripting commands.
Use another program like Avogadro to draw and optimize the geometry and export to read into PyMol later.


(Disclaimer: I'm an Avogadro and Open Babel developer, but there are many programs that could be used.)