New functionality of psikit #Chemoinformatics #RDKit #Psi4

Happy new year! ‘Akemashite Omedetou Gozaimasu’ in Japanese!

I and @kzfm_ -san are developing a library which uses a rdkit & psi4 named psikit. You can find brief introduction the concept following URL.
http://blog.kzfmix.com/entry/1536978824

Today I added new function that can calculate RESP charge of give molecule. RESP charge means ‘Restraints for Deriving
Atomic Charges’. Details can read from this link.

OK, I would like to show how to call new function. It is easy to use but it needs long time for large molecule. I used acetic acid as example. To calculate RESP Charge, user need to run optimization.

############################################
# This code came from RESP branch, not master branch.
############################################

from psikit import Psikit
smi = 'CC(=O)O'
# load smiles and perform optimzation
pk = Psikit()
pk.read_from_smiles(smi)
%time pk.optimize()
>Optimizer: Optimization complete!
>CPU times: user 21.6 s, sys: 920 ms, total: 22.5 s
>Wall time: 6.23 s
>-227.82082530474275

Now I can get HOMO, LUMO and RESP Charge. RESP charges are set as atomic properties of pk.mol object. Let’s check before calculation and after calculation.

# Get Atoms from pk.mol object and check atom properties
atoms = pk.mol.GetAtoms()
print(list(atoms[0].GetPropNames()))
> []
# run the calculation
charges = pk.resp_charge
atoms = pk.mol.GetAtoms()
print(list(atoms[0].GetPropNames()))
> ['EP_C', 'RESP_C']  # EP_C means electro static potential, RESP_C means RESP charge.

Finally draw molecule with RESP charges.

from rdkit.Chem.Draw import MolDrawOptions
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG
view = rdMolDraw2D.MolDraw2DSVG(600,800)
op = view.drawOptions()
for idx, atom in enumerate(pk.mol.GetAtoms()):
    num = float(atom.GetProp('RESP_C'))
    op.atomLabels[idx] = "{}({:,.2f})".format(atom.GetSymbol(), num)
AllChem.Compute2DCoords(pk.mol)
view.DrawMolecule(pk.mol)
view.FinishDrawing()
svg = view.GetDrawingText().replace('svg:', '')
SVG(svg)
acetic acid with RESP charge 6-31G** basis

The results indicates Oxygen atom has the most negative charge and alpha carbon shows more negative than carbonyl carbon.
The library is now under development. We would like to release the package near the feature.

All code of the post can get from the URL below.

https://nbviewer.jupyter.org/github/Mishima-syk/psikit/blob/resp/psikit/example.ipynb

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.