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

Nice post. Thanks!

The transformation matrix is 4×4 in order to allow it to capture rotations and translations.

There’s an explanation here: http://www.cs.cmu.edu/afs/cs/academic/class/15462-s10/www/lec-slides/lec04.pdf

-greg

That makes sense. I understood!!! Thank you for sharing useful information. It is really useful for me. ;)

I’ve come across this here now and working towards a solution. I don’t think cartesian coordinates is the way to do this. I would go back a bit and look into z-matrices.