I was wondering if there were any reliable algorithms out there that could allocate bond orders and formal charges of a given molecular structure if just the atomic coordinates (x,y,z) and atom type were present in a file (C, N, O, H, etc...). The atoms are most likely ones found in organic molecules so complicated heavy metals with multiple charged states are unlikely.

In a step within a molecular modelling workflow, I only have a set of atomic coordinates but no explicit bond or charge information is kept so I need re-derive them as the coordinates have changed (within reasonable extents like bond lengths) relative to the ones from the input structure which has the bond and charge information. So taking the bond and charge information from the input structure probably won't be feasible.

Any references to papers with code would be really helpful for this.

Thank you!

  • $\begingroup$ I would guess any 3D molecular visualization program (e.g. Avogadro, PyMol, Maestro) would have some algorithm for this in order to display molecules in a chemically intuitive way. Avogadro and PyMol both have public repos on Github, so a good place to start might just be searching through the code there. $\endgroup$
    – Tyberius
    Feb 22 at 17:37
  • $\begingroup$ @Joshua I have made an answer. Please check if it answers your question. If it does please click on the tick mark next to it. You also get +2 if you do it. Also acknowledges the effort that goes into answering a question. $\endgroup$ Mar 13 at 16:23

2 Answers 2


This can be achieved by a 2-step process:

Step 1: Convert XYZ to SDF using this open source python library xyz2mol

Step 2: Use the SDF to get the bond order and formal charge information using Python + rdkit.

Below is a fully reproducible example of both steps using Methanol as a sample compound. All of the following can be run in a Google Colab notebook:

!git clone https://github.com/jensengroup/xyz2mol
%cd xyz2mol
!pip install -r requirements.txt
!pip install rdkit

Download the XYZ file of methanol from here.

Convert the XYZ to SDF like this:

!python xyz2mol.py path-to-methanol.xyz -o sdf > path-to-methanol.sdf

Use the SDF for extracting bond order and formal charge:

from rdkit.Chem.rdmolfiles import SDMolSupplier
from rdkit.Chem.rdmolops import GetFormalCharge

suppl = SDMolSupplier('/content/drive/MyDrive/output/methanol.sdf',removeHs=False)
for mol in suppl:
   if mol is not None: 

     for bond in mol.GetBonds():
       print(bond.GetBeginAtom().GetSymbol(),  bond.GetEndAtom().GetSymbol(), bond.GetBeginAtomIdx(), bond.GetEndAtomIdx(), bond.GetBondTypeAsDouble()) 

>> This will output the following:
    C H 3 0 1.0
    H C 1 3 1.0
    H C 2 3 1.0
    C O 3 4 1.0
    O H 4 5 1.0

Here the last column is the bond order.

And to get formal charge :


>> 0

Documentation on the GetBondTypeAsDouble function is available on rdkit here

  • 1
    $\begingroup$ Cool this looks promising. Unfortunately I can't use it directly because this is part of a bigger Rust codebase but I'll try to adapt the approach to the best I can! Thank you:D $\endgroup$ Mar 21 at 10:09

This depends on your workflow purpose. It is usually straightforward to draw chemically sensible bonds between atoms -- this can be done in VMD using Topotools, as one example -- but partial charges are entirely different.

If you simply want a sensible charge distribution for a snapshot somehow, you could try a single structure charge equilibration by loading your structure into LAMMPS and applying a suitable ReaxFF charge equilibration procedure (simply run 0 and write out the data file which will have updated charges).

But you may need transferable charges (every methanol O has the same charge, every acetonitrile N has the same charge, ...) for continuing on to do molecular dynamics (indeed, a good initial charge guess is also very helpful for the above mentioned ReaxFF procedure). In that case you're better off identifying the chemical species in your simulation and then going to either published force fields or something like the Automated Topology Builder for your force field.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .