Calculate atomic charges with psikit #RDKit #psi4

Recently I implemented new function in psikit for atomic charge calculation. Now user can get mulliken charges, RESP charges and Lowdin charges very easily. There is quick example below.

At first import libraries.

import os
import sys
from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
import pandas as pd
import numpy as np
from collections import defaultdict
sys.path.append('/Users/iwatobipen/develop/chemoenv/psikit/psikit')
from psikit import Psikit

I used methanol and imidazole for demo. And psikit can generate 3D structure from SMILES and optimize the conformer with RDKit. Then optimize geometry with psi4.

First example is methanol

pk = Psikit()
methanol = 'CO'
imidazole = 'C1=CN=CN1'
pk.read_from_smiles(methanol)
pk.optimize()
# from literature 6-31G*
# C 0.21 0.04
# O 0.61 0.52
# H 0.13-0.16  0.04
# H 0.39 0.35
resp = pk.resp_charges
mulliken = pk.mulliken_charges
res = defaultdict(list)
for i, atom in enumerate(pk.mol.GetAtoms()):
    res['SYMBOL'].append(atom.GetSymbol())
    res['MULLIKEN'].append(mulliken[i])
    res['RESP'].append(resp[i])
    #print('{:.3f} {:.3f} {}'.format(mulliken[i], resp[i], atom.GetSymbol()))
methanol_df = pd.DataFrame(res)

Oxygen atom shows negative charge. It seems reasonable I think but the values are far from the reference values.

Second example is imidazole.

pk.read_from_smiles(imidazole)
pk.optimize()

resp = pk.resp_charges
mulliken = pk.mulliken_charges
# from literature 6-31G*
#     MULLIKEN RESP
# C 0.02 0.24 
# H 0.15 0.16 
# N 0.55 0.10
# C 0.21 0.06
# H 0.15 0.12
# C 0.03 0.05
# H 0.14 0.11
# N 0.42 0.42
# H 0.33 0.26

#Basis set 6-31G**
res = defaultdict(list)
for i, atom in enumerate(pk.mol.GetAtoms()):
    res['SYMBOL'].append(atom.GetSymbol())
    res['MULLIKEN'].append(mulliken[i])
    res['RESP'].append(resp[i])
    #print('{:.3f} {:.3f} {}'.format(mulliken[i], resp[i], atom.GetSymbol()))
imidazole_df = pd.DataFrame(res)

Two nitrogen atoms shows negative values in both Mulliken/RESP methods.

Reference data came from the URL. I searched mulliken charges of both molecules in google. I found some reports about that but all values are slightly different. I would like to know the accurate value of atomic charges of these molecules. And would like to know how to calculate spin multiplicity with RDKit. Any advices are greatly appreciated. ;)

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.