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. ;)