Relation ship between dihedral deg and atomic charge #psi4 #RDKit #psikit

Recently psikit repository got PR about RESP charge calculation. Thanks for PR. And I have question about the relation ship between compound conformation and partial charge.

Fortunately, psikit already has an example for torsion scan thank @fmkz___ for sharing useful code. The example code is here.

Following code is same as example code linked above. But I calculated not only energy but also RESP charge at each conformer.

Here is an code.

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Draw
from rdkit.Chem import rdMolTransforms
from rdkit.Chem.Draw import IPythonConsole

def mol_with_atom_index(mol):
    for atom in mol.GetAtoms():
    return mol

import py3Dmol

I used 2-Propenoic acid methyl ester as an example.

mol = Chem.MolFromSmiles('C=CC(=O)OC')
hmol = Chem.AddHs(mol)

I generated conformers which have different dihedral angle of bond between carbonyl carbon and carbon of double bond. To generate conformer which has different dihedral angle can generate with SetDihedralDeg method of rdkit. B3LYP/6-31G* is used in following example.

from psikit import Psikit
pk = Psikit()

resp_charges = []
dihedral_energies = []
dihedral_degrees = [i for i in range(0, 360, 10)]

for deg in dihedral_degrees:
    conformer = pk.mol.GetConformer(0)
    rdMolTransforms.SetDihedralDeg(conformer, 0, 1, 2, 3, deg)
    degy ="b3lyp/6-31g*")
    print(deg, degy)
        resp_chage = pk.calc_resp_charges()
        resp_chage = []

OK let’s check results. At first check relative energy of the molecule with different angle.

import numpy as np
dihedral_energies_array = np.array(dihedral_energies)

rerative = (dihedral_energies_array - dihedral_energies_array.min()) * pk.psi4.constants.hartree2kcalmol
%matplotlib inline
import matplotlib.pyplot as plt'ggplot')
plt.plot(dihedral_degrees, rerative)
plt.xlabel('Dihedral Angle degree')
plt.ylabel('Relative energy kcal/mol')
relative energy

The compound is simple alpha-beta unsaturated ester so minimum energy is shown at angle at 0 and 180 degree seems reasonable.

How about resp charge of beta carbon?

chargeofO = []
chargeofC = []
for respc in resp_charges:
        oc = respc[3]
        cc = respc[0]
        oc = None
        cc = None
plt.plot(dihedral_degrees, chargeofC, c='b')
plt.xlabel('Dihedral Angle degree')
plt.ylabel('RESP Charge of Caron_0')

Then I got following figure about carbon which index is zero (beta carbon).

RESP charge of beta carbon at different dihedral angle.

The trend of resp charge of beta-carbon is totally opposite, 0 and 180 deg shows highest charge (90 deg shows most negative value). I think it indicates that electron of conformers which have dihedral angle 0, 180 deg are de-localized. And it is important that partial charge of atoms are strongly depends on molecular conformation. So conformation search and geometry optimization is important task for computational chemistry.

Today’s code can be found from following URL. With nbviewer, 3d compound is visible but gist can’t render 3d molecule.

Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

Published by iwatobipen

I'm medicinal chemist in mid size of pharmaceutical company. I love chemoinfo, cording, organic synthesis, my family.

Leave a Reply

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

You are commenting using your 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.

%d bloggers like this: