An algorithm for implicit bond order and charge allocation for atomic coordinates?

Question

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!

Answer 1

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 :

GetFormalCharge(mol)

>> 0

Documentation on the GetBondTypeAsDouble function is available on rdkit here

Answer 2

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.