On reflection, transforming molecules is tricky

I've just spent quite some time trying to get the mirror-image of a chiral ligand in a metal-ligand complex. It took me a good while to figure out how to do it, so here it is for the record.

The solution uses a combination of Open Babel's API and Numpy. Open Babel can calculate the required transformation matrix given the normal to the mirror plane. To get the normal, we just need three points in the plane, from which we can derive two vectors (it doesn't matter which) which lie in the plane. The cross product of the two vectors is a vector that is orthogonal to both (that's a property of cross products), and thus is a normal to the plane.

It should be possible to do all of the maths with Open Babel, but whatever way the matrix3x3 and vector3 classes are implemented, it doesn't translate into Python (at least not without segfaults). So that explains why I've had to convert to Numpy arrays and do the maths there.

	import pybel
	ob = pybel.ob
	import numpy as np
	def getfrag(mol, startid, ans):
	# Get the disconnected fragment starting at startid
	ans.append(startid)
	for nbr in pybel.ob.OBAtomAtomIter(mol.GetAtom(startid)):
	if nbr.GetIdx() not in ans:
	getfrag(mol, nbr.GetIdx(), ans)
	return ans
	def t(v): # Convert vector3 to np.array
	return np.array([v.GetX(), v.GetY(), v.GetZ()])
	def u(ar): # Convert np.array to vector3
	return ob.vector3(float(ar[0]), float(ar[1]), float(ar[2]))
	def m(mat): # Convert matrix3x3 to np.array
	ans = np.zeros((3,3), "d")
	for i in range(3):
	for j in range(3):
	ans[i,j] = mat.Get(i, j)
	return ans
	if __name__ == "__main__":
	mol = pybel.readfile("out", "CloseToTS.out").next().OBMol
	# Split off the ligand
	bond = mol.GetBond(mol.GetAtom(12), mol.GetAtom(27))
	mol.DeleteBond(bond)
	# Get the indices of atoms in the ligand
	fragment = getfrag(mol, 12, [])
	# Plane composed of Points pqr
	p = t(mol.GetAtom(12).GetVector())
	q = t(mol.GetAtom(50).GetVector())
	r = t(mol.GetAtom(51).GetVector())
	# Normal to plane pqr
	planenormal = u(np.cross(p-q, p-r))
	# Get transformation matrix
	matrix = ob.matrix3x3()
	matrix.PlaneReflection(planenormal)
	mmat = m(matrix)
	# Apply the transformation
	for atom in fragment:
	vec = t(mol.GetAtom(atom).GetVector())
	vec = np.dot(vec, mmat)
	mol.GetAtom(atom).SetVector(u(vec))
	# Add back the ligand bond
	mol.AddBond(12, 27, 1)
	pybel.Molecule(mol).write("sdf", "tmp.sdf", overwrite=True)

view raw mirror.py hosted with ❤ by GitHub