Atomic Charges from OpenBabel

Atomic charges are one of the key parameters that you need to generate for modeling molecular systems. Atomic charges determine the overall electrostatic interaction and are key for modeling intermolecular interactions between protein/nucleic acid/ligand systems.

In force fields, all the parameters are pre-computed except for atomic charges of organic drug-like molecules because of the vast diversity in chemical space that is ever-expanding. Charge models such as CM1A, CM3P, CM1A-BCC, AM1-BCC, RESP, CM5, and CM1A-LBCC use quantum mechanical calculations at semi-empirical or DFT/HF levels of theory for calculating atomic charges, making this step a key bottleneck in molecular simulation setups.

There are charge models such as MMFF94, Gasteiger, and EEM charges, and many more summarized in this J. Cheminf paper that provides quick alternatives for charge generation. These charges can be used in combination with existing force fields such as Amber/GAFF, OPLS-AA force fields to set up systems quickly but their accuracy needs to be investigated. For simple minimizations, these charges are great and are often used in docking programs.

In this blog, I want to share code for generating charges from OpenBabel. This blog is inspired by a really cool Kaggle notebook written by ALEXANDRE SAUVÉ a couple of years ago.

Step 1: Generate 3D molecule from Smiles code

To demonstrate the charge generation from scratch, I will first show how to create a 3d molecule in OpenBabel from smiles code using celecoxib as an example.

import pybel
import openbabel as ob

def Create3DFromSmiles(smiles_code):
    mol = ob.OBMol()
    conv = ob.OBConversion()
    conv.SetInFormat('smi')
    conv.ReadString(mol, smiles_code)
    pbmol = pybel.Molecule(mol)
    pbmol.addh()
    pbmol.make3D(forcefield="gaff",steps=100)
    pbmol.localopt()
    pbmol.write('pdb','molecule.pdb',overwrite=True)
    return pbmol.OBMol

This block takes smiles code as input and generates molecule.pdb with 3d structure.

mol = Create3DFromSmiles('O=S(=O)(c3ccc(n1nc(cc1c2ccc(cc2)C)C(F)(F)F)cc3)N')

Running it successfully generates a structure that looks like this in ChimeraX

ChimeraX Stick Representation of Celecoxib

Step 2: Calculate charges using the EEM charge model

ob_charge_model = ob.OBChargeModel.FindType("eem2015bn")
ob_charge_model.ComputeCharges(mol)
charges = ob_charge_model.GetPartialCharges()

The block above generates an array of EEM charges. Now write them as b-factors and visualize them using ChimeraX

from biopandas.pdb import PandasPdb
ppdb = PandasPdb()
ppdb.read_pdb('./molecule.pdb')
ppdb.df['HETATM']['b_factor'] = charges
ppdb.to_pdb(path='./molecule_wQs.pdb',records=None,gz=False,append_newline=True)

Step 3: Color atoms by charge value in ChimeraX

Open molecule_wQs.pdb in ChimeraX and type the following commands to generate the figure below. You will need to know the minimum and maximum values of charges generated for the color bar.

color by bfactor #1 range -1,2
key blue:-1 white: red:2
save image_wQs.png format png supersample 4 transparentBackground true

Celecoxib with atoms colored by EEM2015bn charges calculated from OpenBabel

There are a variety of other charge models that you can generate from openbabel such as the ones below.

charge_models = [ "eem", "mmff94", "gasteiger", "qeq", "qtpie", 
                  "eem2015ha", "eem2015hm", "eem2015hn", 
                  "eem2015ba", "eem2015bm", "eem2015bn" ]

The notebook I referenced at the beginning of the blog nicely summarizes which ones are the most similar to Mulliken charges. eem&eem2015bnhave the highest R2 value with Mulliken charges and hence I used them in this blog for demonstration. Feel free to play around with other methods. Hope this is useful for you.

References

[embed]V7 Estimation of Mulliken Charges with Open Babel Explore and run machine learning code with Kaggle Notebooks | Using data from Predicting Molecular Propertieswww.kaggle.com

[embed]Command: color, rainbow Depending on the type of coloring, color may apply to different sets of items: Simple Coloring - atomic models,volume…www.cgl.ucsf.edu

[embed]Command: key The key command creates a legend or key to show the value-color correspondences in a figure, such as from: It is the…www.cgl.ucsf.edu

[embed]High-quality and universal empirical atomic charges for chemoinformatics applications - Journal of… Background Partial atomic charges describe the distribution of electron density in a molecule and therefore provide…jcheminf.biomedcentral.com

Shameless self-plug

[embed]*1.14CM1A-LBCC: Localized Bond-Charge Corrected CM1A Charges for Condensed-Phase Simulations* The quality of the 1.14CM1A and 1.20CM5 charge models was evaluated for calculations of free energies of hydration…*pubs.acs.org

[embed]LigParGen Edit descriptionzarbi.chem.yale.edu