Rotate molecule and visualize it #RDKit #Chemoinformatics

Some days ago, I found an interesting question in rdkit mailing list. The question is how to rotate molecule around an axis. I do not have any idea to do it. But RDKit has a function to do it. I read API and try to do it.

rdMolTransforms.TransformConformer function can rotate molecue with transform matrix. So, I think I can do it if rotate matrix is defined.

Rotation matrix is a simple. Like below.

def rot_ar_x(radi):
    return  np.array([[1, 0, 0, 0],
                      [0, np.cos(radi), -np.sin(radi), 0],
                      [0, np.sin(radi), np.cos(radi), 0],
                     [0, 0, 0, 1]], dtype=np.double)

def rot_ar_y(radi):
    return  np.array([[np.cos(radi), 0, np.sin(radi), 0],
                      [0, 1, 0, 0],
                      [-np.sin(radi), 0, np.cos(radi), 0],
                     [0, 0, 0, 1]], dtype=np.double)

def rot_ar_z(radi):
    return  np.array([[np.cos(radi), -np.sin(radi), 0, 0],
                      [np.sin(radi), np.cos(radi), 0, 0],
                      [0, 0, 1, 0],
                     [0, 0, 0, 1]], dtype=np.double)

Let’s write code!

import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import rdMolTransforms
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit import rdBase
from IPython import display
import copy
import py3Dmol

def rot_ar_x(radi):
    return  np.array([[1, 0, 0, 0],
                      [0, np.cos(radi), -np.sin(radi), 0],
                      [0, np.sin(radi), np.cos(radi), 0],
                     [0, 0, 0, 1]], dtype=np.double)

def rot_ar_y(radi):
    return  np.array([[np.cos(radi), 0, np.sin(radi), 0],
                      [0, 1, 0, 0],
                      [-np.sin(radi), 0, np.cos(radi), 0],
                     [0, 0, 0, 1]], dtype=np.double)

def rot_ar_z(radi):
    return  np.array([[np.cos(radi), -np.sin(radi), 0, 0],
                      [np.sin(radi), np.cos(radi), 0, 0],
                      [0, 0, 1, 0],
                     [0, 0, 0, 1]], dtype=np.double)
tforms = {0: rot_ar_x, 1: rot_ar_y, 2: rot_ar_z}

I used simple molecule for test.

mol = Chem.AddHs(Chem.MolFromSmiles("COC1CN(C)CC1"))
AllChem.EmbedMolecule(mol)

Then I defined draw function. The function borrowed from rdkit-blog. ;)

# ref  http://rdkit.blogspot.com/2016/07/using-ipywidgets-and-py3dmol-to-browse.html
def drawit2(m,p,confId=-1):
    mb = Chem.MolToMolBlock(m,confId=confId)
    p.addModel(mb,'sdf')
    p.setStyle({'stick':{}})
    p.setBackgroundColor('0xeeeeee')
    p.zoomTo()

Now ready. Let’s rotate molecule around x, y, z axis.

# around X axis
p = py3Dmol.view(width=400, height=400)
for i in range(11):
    rdMolTransforms.TransformConformer(mol.GetConformer(0), tforms[0](2*np.pi/10))
    drawit2(mol,p)
p.show()

Looks good.

# Before align
v = PyMol.MolViewer()
v.DeleteAll()
process = subprocess.Popen(['python', showfeatpath, '--writeFeats','./before_align_cdk2.sdf'], stdout=subprocess.PIPE)
stdout = process.communicate()[0]
png=v.GetPNG()
display.display(png)

Rotate Z axis last!

# around Y axis
p = py3Dmol.view(width=400, height=400)
for i in range(11):
    rdMolTransforms.TransformConformer(mol.GetConformer(0), tforms[1](2*np.pi/10))
    drawit2(mol,p)
p.show()

Works fine. I do not understand why transformation matrix is 4 x 4. I wonder if you have any comments about it.

The code is uploaded to my repo and you can check from following URL.
https://nbviewer.jupyter.org/github/iwatobipen/playground/blob/master/rotation_mol2.ipynbh

Advertisement