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
AllChem.SetPreferCoordGen(True)

def mol_with_atom_index(mol):
    for atom in mol.GetAtoms():
        atom.SetAtomMapNum(atom.GetIdx())
    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)]
pk.read_from_smiles('C=CC(=O)OC')
pk.optimize(basis_sets="scf/sto-3g")
pk.optimize(basis_sets="b3lyp/6-31G*")

for deg in dihedral_degrees:
    conformer = pk.mol.GetConformer(0)
    rdMolTransforms.SetDihedralDeg(conformer, 0, 1, 2, 3, deg)
    degy = pk.energy(basis_sets="b3lyp/6-31g*")
    print(deg, degy)
    try:
        resp_chage = pk.calc_resp_charges()
    except:
        resp_chage = []
    dihedral_energies.append(degy)
    resp_charges.append(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
plt.style.use('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:
    try:
        oc = respc[3]
        cc = respc[0]
    except:
        oc = None
        cc = None
    chargeofO.append(oc)
    chargeofC.append(cc)
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.

https://nbviewer.jupyter.org/gist/iwatobipen/3e7bc08a0872c44a251c5c42cada0011

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw torsion_resp.ipynb hosted with ❤ by GitHub

Advertisement